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EXTENSION  OF  PENCIL  OF  FUNCTIONS  METHOD  TO  REVERSE-TIME 
PROCESSING  WITH  FIRST-ORDER  DIGITAL  FILTERS 


I.  INTRODUCTION 

Signal  representation  and  approximation  is  basic  to  (a)  time-domain 
extraction  of  singularities  of  a  scatterer's  field  pattern.  It  is  also 
useful  in  (b)  bandwidth  compression  of  signals,  and  (c)  time-domain  measure¬ 
ment  and  testing  of  networks /channels.  This  report  discusses  a  unified 
approach  to  representing  or  approximating  a  given  empirical  signal  by  sum 
of  exponentials,  i.e.,  for  finding  the  right  hand  side  of 

n  A.t  n  W. 

h  ,(t)  =  h(t)  -  Z  W,  e  1  —  H(s)  -  Z  - 

d  i-1  1  i-1  8  "  Ai 

.  ^n-2s°~^+*  **  8o 

sn  +  a  .+....  +  a 

n-l  o 

or,  equivalently,  the  right  hand  side  of  the  sampled  version 

n  .  n  R 

h.(k)  =  h(k)  -  Z  R,(r.r  ~  H(z)  -  Z  - ^-r- 

d  i-1  1  1  i-1  (l-z±z  X) 

.  ,  ,  -1,  ,  .  -n+1 

b  +  b.z  +  ...+  b  .z 
o  1  n-l _ 

1  +  a. z^+...+azn 
l  n 

The  poles  (or  in  z-domain)  are  either  real,  or  they  occur  in  complex 
conjugate  pairs. 

In  the  method  described  here,  the  data  signal  is  processed  in  reverse- 
tine  by  a  cascade  of  first  order  digital  filters  —  each  y(z)  -  l/(l-qz 
to  yield  a  family  of  information  signals.  The  Gram  matrix  Q  of  these  infor¬ 
mation  signals  is  shown  to  contain  the  essetltial  information  on  the  denominator 
parameters  of  H(z).  Specifically,  it  is  shown  that  A(z)  is  determined  as 
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A(z)  -  z~n  [  E  /D.  (z-l)n+1“iJ//D, 
i-1  1  1 

where  are  the  diagonal  cofactors  of  the  matrix  G«  The  numerator  para¬ 
meters  are  then  determined  using  a  least-squares  fit,  i.e.,  b_  *  -  P”^v» 
where  P  and  v  are  defined  in  the  paper. 

The  entire  procedure  is  thus  noniterative  and  computationally  efficient. 
It  is  a  further  generalization  of  the  method  developed  in  [4  ] .  Examples 
presented  demonstrate  (i)  noiseworthiness  in  the  representation  problem 
when  data  are  corrupted  by  noise  and  (ii)  the  effectiveness  of  the  method 
in  the  approximation  problem.  Comparison  of  the  method  with  the  maximum 
entropy  method  (or  linear  predictor)  and  the  Prony  method  is  also  included 
in  the  report. 

The  structure  of  the  report  is  as  follows.  The  fundamental  results 
relating  to  the  signal  representation/approximation  problem,  through  the 
use  of  reverse-time  processing  by  first-order  digital  filters,  are  obtained 
in  Chapter  II.  The  important  question  of  parallelism  and  orthogonality  of 
the  Information  signals  is  explored  in  Section  III.  Chapter  IV  gives  the 
input  description  and  user  instructions  for  the  program  POF-FILTER  which 
Implements  the  method  of  Chapter  II.  Application  examples  as  well  as  com¬ 
parison  with  other  methods,  are  given  in  Section  V.  Appendices  A  gives 
listing  and  description  of  program  POF-FILTER  and  its  routines.  Appendices 
B  and  C  contain  the  computer  outputs  relating  to  Example  1  and  Example  2, 
respectively,  of  Section  V. 


2 


SECTION  II 

FIRST-ORDER  FILTER  BASED  PENCIL-OF- FUNCTIONS 
METHOD  FOR  MODELING  IMPULSE  RESPONSES 


We  shall  be  interested  in  modeling  the  impulse  resonse  [l]-[3] 


y(k)  -  I  R?  (Zn) 
£*  1 


.  b.  z"1  +  ...  +  b  z“n 

Mil  .  _ i _ S - 

A(z)  1  +  a,  z_1+  . . .  +  a  z“n 
1  n 


(1) 


from  its  numerical  data.  Suppose  a  suitable  K  has  been  selected  such  that 
y(k)  ■  0  for  k  >  K  (so  that  use  of  the  upper  limit  K  instead  of  00  on  summa¬ 
tions  may  be  permitted).  We  define  the  reverse-time  first-order  filtered 
signals  as 


y^k)  *  y(fc) 


y2<k)  -  qy2(k+l)  +  y^k) 


(2) 


yN(k)  -  qyN(k+l)  +  yn(k) 

where  N-n+1,  and  y^(K+l)-0  for  i-l,2,..,N.  Further,  0  <  q  <  1. 

This  family  of  signals,  which  we  shall  call  information  signals, 
exhibits  the  interesting  property  stated  below. 

Lemma  1 


,  (k)  -  l 


R„ 


(3) 


1+1  Jt-1  (1-q  z^)1 

Proof;  We  prove  this  by  induction.  For  i“0  the  statement  is  trivially 
true  since  it  is  identical  to  (1)  for  this  case.  Assuming  it  to  be  true 
for  i-1,  let  us  proceed  to  prove  it  is  true  for  i. 

From  (2) 


y1+1(k)  -  qyi+1(k+l)  +  yt(k) 


(A) 
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which  is  readily  shown  to  be  equivalent  to 


Wk) 


£  yt(v) 

v»k  1 


n 

Z 


£-1  (1-q  z.)i_1  V-k 


v  v-k.  .v 

2  q  (z0) 


(from  induction 
hypothesis) 


n 

-  Z 


R, 


k  _  v-k  .  Nv-k 

—  *.  Z  q  (zJ 

1  1  v-k  * 


£-1  (1-q  z^) 

The  result  of  equation  (3)  follows  immediately  by  observing  that 

oa 

1 


r  v— k  ,  x V-k 

2  q  (zff) 

v-k 


U-q  z£) 

This  Lemma  leads  us  to  the  crucial  observation  stated  next. 
Lemma  2 

The  set 

(qzn,-l)y2  +  yr  (qzm-i)y3  +  y2,  ...  ,  (qV1)yN  +  yn 


(5) 


is  linearly  dependent  for  m-l,2,...,n  where  z  are  the  poles  of  the  right 
•  in 

hand  side  of  (1). 

Note:  We  have  used  the  notation  y^  to  denote  the  sequence  (y^(k)}, 
k-0 ,1,2 . 

Proof:  In  view  of  (3)  we  find 


(qv1)yi+i(k)  +  yi(k)  “  } 


n  q(vzA,  xk 

- r-  (z,) 

A- 1  (1-q  z,) 

Wm  * 


(6) 


for  i-l,2,...,n.  Clearly,  the  sequences  Cqz  -Dy^j+y^ ,  i-1,...,  n 
each  contain  only  n-1  modes  (z^)  ,  Jlfm  hence  are  linearly  dependent. 
We  can  now  apply  the  pencil-of-functions  theorem  of  reference  [4] 


to  obtain  the  central  theoretical  result  of  this  section. 


We  will  call  (see  (1))  the  poles  of  the  impulse  response,  the 

I- 

corresponding  residues,  and  p^-  { ( z g )  }  the  associated  modes.  Note  that 
the  poles  occur  in  conjugate  pairs  whenever  complex,  as  do  the  residues, 
since  y  is  real. 

Define  the  N  x  N  dimensional  Gram  matrix  (recall,  N-n+1)  [6] 


F  - 


<yl,yl>  *  '  *  <yl,yN> 


•  •  •  <yN'V 


K 


<y.,y,>  -  l  y, (k)y  (k)  (7a) 

1  3  k-1  1  1 


or,  equivalently, 

K 

F  -  Z  f(k)fT(k) 
k-1 

where  (k)  -  [y^k)  y2(k)  ....  yH(k)]. 


(7b) 


Theorem  1 


The  poles  ol  the  Impulse  response  y(fc)  must  satisfy  the  equation 

-  0  (8) 


N  N-i 

Z  /D  (qz-1) 
i-1  1 


where  are  the  diagonal  cofactors  of  the  Gram  matrix  F  (defined  in  (7) 
Proof:  The  theorem  follows  inmediately  upon  application  of  the  pencil 
of  functions  theorem  (reference  [4])  to  the  set  (5). 

Note  that  the  denominator  of  transform  of  the  impulse  response 
is  given  by 

A(z)  -  D^172  (qz)”n  Z  A>  (qz-l)N*i  (9) 

i«l 

This  follows  from  (8)  by  dividing  through  by  zn  and  by  normalizing  the 
coefficients  so  that  the  leading  coefficient  becomes  unity. 
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Determination  of  the  denominator  polynomial  A(z)  completes  the  first 


of  the  two  steps  in  the  decoupled  pencil-of-f unctions  method.  The  next 
step,  finding  the  numerator  coefficients  in  B(z),  or  equivalently  finding 
the  residues  R^,  can  be  accomplished  in  two  alternative  ways. 

The  first  method  consists  in  solving  for  the  residues  from  the 
equation  ^ 


1 

1  .  . 

1 

D 

l-qzj^ 

i-qz2 

1-qz 

n 

R1 

1 

1  .  . 

i 

1  j 

(l-qz1)1 2 * * * 6 

• 

(l-qz2)2 

•  «  • 

(1-q*  )2  ! 
n 

i 

* 

R 

R2 

• 

• 

1 

•  •  • 

1 

i 

. 

I 

i 

1 

• 

(1-qZj) 

(l-qz2) 

/i  \N 

(l-qzn> 

K 

“ 

y2(o) 


y3(0) 


(10) 


yN(°) , 


This  equation  follows  from  (3)  upon  setting  k«0  and  letting  i  range  from 
1  to  n.  Clearly,  the  use  of  this  equation  requires  that  the  poles  z  be 
determined  from  the  denominator  A(z)  by  use  of  a  root-finding  routine. 

This  requirement  is  of  no  consequence  if  the  final  answers  are  needed  in 
the  s-domaln,  for  conversion  to  s-domain  involves  finding  the  roots  of 
A(z)  anyway. 

The  alternative  approach  is  to  find  the  optimum  least-squares 
numerator  coefficients  (given  the  denominator  of  (9))  through  the  equation 


1  The  equation  corresponding  to  i«0  is  ignored  in  formulating  (10)  because 

of  the  relatively  poor  signal/noise  statistics  of  y^(0)  (compared  to, 

say,  y.(0))  when  the  impulse  response  data  contain  additive  wideband 

noise.  This  statement  assumes  that  the  bandwidth  of  y(k)  is  much 

smaller  than  that  of  the  additive  noise. 
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A 


<wr 

V 

<wr 

w2> 

•  •  •  <w ^ | 

w  > 
n 

bl 

A 

* 

*-» 

V 

<W2* 

wl> 

<W2* 

V 

•  •  •  ^2 1 

w  > 
n 

b2 

<y.  »2> 

• 

• 

•  •  •  • 

i 

• 

• 

(id 

JV 

wl> 

<Wn* 

V 

•  •  •  • 

.  .  •  <w  , 
n 

! 

w  >J 

n  J 

• 

Lv 

1  • 

L<y.  V- 

where  w^  denotes  the  impulse  response  of  z  */A(z).  Note  that  wi(k)»w(k-i) 
where  w(k)  is  the  impulse  response  (i.e. ,  inverse  z-transform)  of  1/A(z). 
All  inner  products  are  summed  from  k*0  to  K. 

Discussion 

Equations  (9)  and  (11)  have  been  Implemented  in  a  computer  program 
"POF-FILTER"  written  in  FORTRAN  IV.  The  program  is  presented  in  Section 
III. 

The  idea  of  reverse-time  Integration  was  proposed  by  Carr  in  [7] 
and  Jain  in  [8].  Here,  we  have  generalized  the  concept  of  reverse-time 
processing  to  the  case  of  first-order  filter  processing.  Note  that  the 
first  order  filter  l/(l-qz  ),  used  above,  encompasses  Integration;  just 
let  q-1. 

It  should  be  borne  in  mind  that  the  approach  developed  above  is 
applicable  to  impulse  responses  only.  With  some  effort,  it  may  be  modified 
for  use  to  step  responses  and  square-pulse  responses.  For  more  general 
Inputs,  however,  one  must  use  the  coupled  approach  dlcussed  in  [9]  (which 
of  course  involves  greatly  increased  computations,  e.g.,  the  Gram  matrix 
involved  is  of  an  order  twice  as  high  as  in  the  decoupled  procedure) . 

The  decoupled  approach  can  be  used  only  if  the  data  is  of  sufficient 

length  K  such  that  y(k>K)-0  for  all  practical  purposes. 
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The  reverse-cime  processing  of  y(k)  through  the  first  order  filters 
can  he  interpreted  as  forward-time  processing  of  the  signal  h(k)*y(K-k) 
through  the  same  filters. 

The  use  of  square-roots  of  cofactors  of  the  Gram  matrix  is  analogous 
to  the  use  of  square-root  factorization  [ 10]  of  the  Gram  matrix.  Attendant 
advantages  are  therefore  expected  to  be  realized.  A  more  detailed  analysis 
of  this  connection  will  be  discussed  elsewhere. 

The  transfer  function  of  the  first  order  filter  used  in  equation  (2) 
is  u(z)»  l/(l-qz-*).  Instead,  we  could  use  filters  with  transfer  funnction 
y^(z)*  (l-q)/(l-qz-*) ;  these  filters  have  a  d.c.  gain  equal  to  unity  and 
the  ratio  of  the  output  power  to  the  input  power  is  a  direct  measure  of 
the  extent  of  the  rejection  of  higher  frequencies.  Equation  (9)  remains 
valid  even  when  these  unity  d.c.  gain  filters  are  used. 

Note  that  the  first-order  filtering  (in  reverse  time)  is  achieved 
in  (2)  recursively,  without  the  need  to  carry  out  discrete  convolution. 
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SECTION  III 

PARALLELISM  AMD  ORTHOGONALITY  OF  INFORMATION  SIGNALS 


Here,  we  consider  Che  important  matter  of  parallelism  and  orthogonality 
of  the  information  signals.  In  the  last  section  we  processed  the  signal 
h(k)«  y(K-k)  through  first  order  filters  P^CzWl-qJ/Cl-qz”*).  This  is 
depicted  in  Fig.  1.  Note  that  forward  time  processing  of  h(k)  is  equivalent 
to  reverse-time  processing  of  y(k).  Because  of  familiarity  with  forward¬ 
time  processing  we  will  carry  out  the  discussion  below  in  terms  of  the 
signal  h(k).  .  Also  note  that  h^(k)»y^(K-k) ,  i«l,2,..,N.  Finally,  we 
remark  that  the  Gram  matrix  and  the  properties  of  parallelism/orthogonality 

of  the  two  families  of  signals  h^.hj . h^  and  yj/Fj*  *  ’  *  ,yN  are  Identical. 

Indeed,  <1^,  h±>  -  <yit  y±>  and  <1^,  h^>  »  <yi>  y^>  for  i,j-l,2 . N. 


h(k)-h  (k)  h,  •  h,  1  K 

- Ak'y^(z)  —  •  >  ...  -UjCz)-  > 


Fig.  1.  Generation  of  information  signals 
by  use  of  first-order  filters. 

The  magnitude  (vs.  frequency)  characteristic  of  the  first  order  filter 
y. (z)  is  shown  in  Fig.  2.  The  cutoff  frequency  (-  3dB  point)  to  is  related 

L  C 

to  the  parameter  q  as  [5] 

«c  ■  -r  >*<t>  <l2> 

where  A  is  the  sampling  interval.  Defining  SI  as  the  normalized  frequency 

(i.e.,  the  ratio  co/to  where  to  is  the  sampling  frequency  in  radian/s) 
s  s 

we  can  express  the  normalized  cutoff  frequency  as 

nc“'HLn(1/q>  (l3) 
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Fig.  3  Normalized  cutoff  frequency  vs.  paroeter  q 


The  normalized  cutoff  frequency  flc  is  plotted  against  the  filter 
parameter  q  in  Fig.  3. 

To  analyze  the  parallelism/ '"thogonality  (PO)  properties  of  the 
information  signals  let  us  define  the  connection  filters 


Attenuation 


so  that  we  may  write 


Hi+l(z)  “  Mi(z)  Hl(z) 


Note  that  M^(z)  is  an  i-pole  filter  with  a  pole  of  multiplicity  i  at  z*q. 
The  magnitude  characteristic  of  the  filters  Mq,...,M^  are  shown  in  Fig.  4 
for  q»0.8. 


.rnwunmaiioui 

liwiVjfiaVHBI 

anvaamm 


0.025 


-20  i - 1 — i - 1— j  i  '  1  i  t  i  i 1  \^vrrr  \  ■  \  ■  i  i.:i  n  o.oi 

0.  lu  l.Ow  *  10. 0u> 

c  c  c 

Fig.  4.  Magnitude  vs.  Freq.  of  connection  filters 
Ue  now  present  an  approximate,  and  somewhat  heuristic,  analysis  of 

the  PO  properties  of  the  family  of  information  signals. 

Approximations  - 

He  approximate  the  magnitude  characteristics  of  the  connection 
filters  as  shown  in  Fig.  5a. 

The  phase  properties  of  will  be  ignored  (assumed  identically  zero) 
It  is  then  possible  to  write 


M,  (e^*1^)  ■  L,  (oi)  +  ....+  L  (u) 
li  n 


where  L^(u>)  are  the  bandpass  filter  characteristics  shown  in  Fig.  5b. 
Further,  (15)  and  (16)  together  yield 
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Power  level 


-10  dB 


U4  “3  W2  “l 


-10  dB 


L0(w) 


Fig.  5.  Idealized  connection  filters  and  BP  constituents. 
Hi+i(ejaiA)  -  (Li (u>)  +  ....  +  Ln(w»  Hl(eJ“&) 


(17a) 


where 


■  ♦ .  (u)  +  ....+♦  (w) 

i  n 


^(w)  -  L±(co)  HjCe^) 


(17b) 


Clearly,  ^(oi),  i-0,l,...,n  are  orthogonal.  Indeed, 

V2  ^  rV2 

J  *1(u)  *j*(»)  du  -  } Hx (oj)  | 2  L^otfLj^w) 


-“s/2 


-w8/2 

-  0  for  i?*j 


because  L^(u)»0  wherever  (u)iH)  and  vice-versa.  Let  us  state  this  in 
the  tine  domain  as 


«J>t,  ^  > 


aii  5ij 


—  * 

where  is  the  Kronecker  delta  and  {^(k)}  -J  *i(u). 

We  can  summarize  the  above  discussion  as  follows. 

Theorem  2 

The  information  signals  are  approximately  tri-orthogonal. 
Proof:  From  (18)  and  (20)  we  have 


1 


1 


V 

h2 

h3 

• 

e 

_v 

1 

0 


.  .  1 
.  .  1 

•  .  1 

•  .  1 

•  •  • 

0  0  1 


u'nJ 


(21) 


where  . »  $n  are  a  set  of  orthogonal  signals.  The  approximation 

sign  arises  because  of  the  assumptions  made  earlier. 

Now  reversing  the  time-indices  about  k«K,  we  have 

l  1  l  .  .  l“j 

0  1  1  .  .  1 

0  0  1  .  .  1 

•  •  •  ♦  •  1 

•  •  •  •  •  • 

o  o  o  ooiJL<p 


Ly3. 


where,  by  definition  ^(k)  -  ^(K-k).  Since 


"  <*i,  ipj> 


•  o  6 
1  1J 


the  assertion  of  the  theorem  Is  proved. 


(22) 


(23) 


Note  that  the  weights  Q  can  be  calculated  (In  view  of  (18))  as 

«.n 

°i  "  /  lLi<aj)  Hl<*Ja^)  1 


-w8/2 

Else,  they  may  be  calculated  In  the  time  domain  via  (20). 
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Note  that  (22)  may  be  writen  more  compactly  as 


I  ■  (25) 

where  £-  [jv  y2,  ...  ,yN  3  T,  [aQ,  . aj  T  and  U  is 

the  upper  triangular  matrix  of  l's.  Secondly,  we  observe  that  the  Gram 
matrix  of  the  Information  signals  can  be  written  as 

F  -  <it  ZT>  (26) 

T 

where  the  inner-product  is  taken  over  each  term  of  the  matrix  .  In 
terms  of  the  approximate  tri-orthogonal  characterization  in  (22) ,  or  (25) , 
we  have 

F  s  0  <J,  iT>  B1 

-  U  E  UT  (26) 


where  E  is  the  n+1  dimensional  diagonal  matrix  diag{an,o, ,  ...  ,a  }. 

0  1  n 

The  lemma  below  summarizes  the  PO  properties  of  the  information 
signals . 

Lemma  3 

The  correlation  coefficient  of  a  pair  of  information  signals 
y±  and  yj ,  j  >  i  is 


3  -  /qj-l+  •"  +  °n 

«  "/  Vl+  +qn  ’ 


i,j«l, . . . ,n+l 


(27) 


Proof  -  The  relation  follows  readily  from  the  approximate  tri-orthogonal 
characterization  of  the  information  signals.  Using  (22)  and  (23),  we  have 

Vl+  -  +gn _ 


Pi^  /ot_]+  . ..+  Oq  /  0^+  ...+ 


which  is  the  same  as  (27). 
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Discussion 

In  any  system  identification  technique  it  is  desirable  to  have  basis 

functions  that  differ  significantly  from  each  other.  Ideally,  they  should 

be  orthogonal.  By  varying  the  parameter  q  between  0  and  1,  it  is  possible 

to  generate  sets  of  basis  functions  that  vary  between  a  set  whose  elements 

are  almost  identical  to  approximately  a  set  whose  elements  are  orthogonal. 

If  q  -  0  (recall  that  q  is  the  z-domain  pole  of  the  first  order 

filters  y^(z)),  then  the  connection  filters  all  have  a  nearly  all-pass 

characteristic,  and  as  a  result  0.  -  0  for  i  n.  Thus  0  dominates  all 

1  n 

other  and  we  have 


This  is  undesirable;  however,  such  is  indeed  the  case  when  unit-delays 
are  employed  such  as  in  methods  like  Prony,  Linear  Predictive  analysis, 
Maximum-Entropy  method,  etc.  [l]-[2].  The  strong  correlation  between  the 
basis  functions  leads  to  numerical  difficulties. 

On  the  other,  extreme,  q-1  results  in  pure  (digital)  integration  of 
y(k) ,  in  reverse-time,  for  generation  of  the  information  signals.  In  this 
case  the  analysis  of  PO  properties  developed  above  is  not  applicable. 
Therefore,  consider  q  close  to  one  (from  the  left;  e.g.,  q  ■  0.99.  Now 
the  cutoff  frequencies  w^,  u^,  ...,  wn  become  crowded  near  the  zero  frequency. 
Hence  0^  =  0  for  i  +  0  and  Oq  dominates  all  other  c^.  Therefore,  .... 
yN  have  very  small  energy  (thus  very  small  fraction  of  signal  informatipn) , 
and  they  are  nearly  orthogonal  to  y^. 

Intermediate  values  of  q  lead  to  other  useful  sets  of  basis  functions. 
Example 

Consider  that  the  signal  y(k)  has  flat  low-pass  spectrum  from  0,  •  0 
to  0.05  with  amplitude  level  10.  Then  q  ■  0.8  and  n  *  4  lead  to  the  values 

oQ  -  1.45,  ol  -  1.28,  o2  -  0.46,  o3  -  0.31,  c ^  -  1.5 

from  which  the  correlation  coefficients  can  be  computed  readily.  For 
example  p,,3  •  0.893,  and  ■  0.548. 
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SECTION  IV 


PROGRAM  DESCRIPTION 

The  program  POP-FILTER  is  a  high  accuracy  FORTRAN  IV  program  designed 
Co  implement  Che  decoupled  pencil-of-functions  method.  Specifically,  it 
models  Che  impulse  response  of  a  finite-order  linear  system  by  processing 
the  given  data  through  a  cascade  of  first-order  linear  filters  in  reverse¬ 
time.  Some  of  the  features  of  the  program  are  stated  below. 

*  Decoupled  denominator  and  numerator  determination.  This  permits 
fairly  high  order  models  to  be  determined,  since  the  order  of  the 
Gram  matrix  is  n+1  where  n  is  the  model  order  (or  n+2  when  data 
bias  is  also  to  be  estimated.  In  contrast,  the  coupled  procedure 
requires  the  use  of  a  2n+2  order  Gram  matrix  (or  2n+3  when  bias  is 
also  to  be  estimated) 

*  First-order  filter,  rather  than  pure  integrators,  are  used  for 
generating  the  family  of  Information  signals.  This  results  in  a 
nearly  tri-orthogonal  set  of  signals  which  result  in  a  better 
conditioned  Gram  matrix  than  with  the  use  of  pure  Integrators. 

*  Bias  correction  option.  Data  bias  can  be  estimated  and  thereby 
more  accurate  estimates  for  the  model  transfer  function  parameters 
obtained. 

*  Noise  correction  option.  A  preliminary  routine  for  estimation  of 
noise  effect, and  correction  thereof,  has  been  included.  Theoretical 
work  and  testing/improvement  of  this  routine  remains  to  be  done. 

*  Direct  transmission  option.  Structural  correctness  of  the  model 
in  either  the  presence  or  absence  of  the  direct  transmission 
term  is  preserved  by  exercising  this  option. 

*  Results  of  modeling  are  obtained  in  the  z-domain  and  on  an  optional 


i-  '  mm 


basis  in  Che  s-domain  also. 

*  Simulation  option.  The  data  modeled  may  be  laboratory  test  data,  or 

one  may  generate  an  impulse  response  within  the  program  by  specifying 

n  k 

it  either  in  the  form  I  R0  Cz„l  or  by  the  z-domain  transfer 

,k-l  *  * 

function  H(z)"(b„+b,z  +  . .+b  z  ;/(l+a,z  +  ..  +a  z  n  ). 

0  1  n  l  n 

In  the  simulation  mode,  a  desired  amount  of  bias  and/or  additive  noise 

/ 

may  be  Incorporated  for  test  purposes. 

In  the  laboratory-data  case,  a  preliminary  bias  removal  and  a  data  scaling 

procedure  (Co  maximize  the  effectiveness  of  the  algorlthm)has  been  incorporated. 

*  The  routine  which  finds  the  cofactors  of  the  Gram  matrix  of  the 
Information  signals  has  been  optimized  by  incorporating  a  scaling 
and  corresponding  descaling  stages. 

*  For  comparison  purposes  the  program  provides  the  option  to  use  two 
other  modeling  techniques,  the  linear  predictor  and  the  prony  method. 

The  latter  can  be  the  classical  Frony,  which  uses  2n+2  data  points, 
or  the  least-squares  Prony. 

The  provision  of  these  methods  within  this  single  program  is  an 
essential  step  coward  the  evaluation  of  the  pencil-of-functions  method 
against  other  benchmark  methods. 

The  input  data  cards  on  the  subsequent  pages  give  a  description  of  all 
input  variables,  and  in  so  doing  provide  an  understanding  of  the  program 
use. 
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INPUT  DATA  CARDS 


CARD  If  1  The  first  card  is  a  title  card. 

The  first  card  is  a  title  card.  Columns  1  through  70  are 
available  for  alphanumeric  title. 

This  title  is  reproduced  in  the  output 

Another  title  card  (columns  1-70) .  Not  reproduced  in  output 
CARD  If  3  Another  title  card  (columns  1-70).  Not  reproduced  in  output 


CARD  If  1 

CARD  If  2 


CARD  If  4 


First  option  card 


Variable  Name 
(Format) 


Description 


Preferred/default 
Columns  Value  _■  .. 


NPT  (14)  Number  of  signal  points  used  in  1-4 

analysis /modeling 

IXX,  IYY  (212)  Blank  (unused)  5-8 

N  (12)  Model  order  9-10 

ISIM  (12)  Simulation  mode  option  11-12 


ISIM  «  -1  Real  (laboratory)  data 
in  integer  format  (1015) 

■  0  Real  (laboratory)  data 

in  real  format  (10F8.6) 

■  1  simulation  from  digital 

transfer  function 
*  2  simulation  from  sums  of 

exponential*sinusoid  form 
(specified  in  continuous-time 
domain) 


NCOMP  (12)  Number  of  terms  in  the  sums  of  13-14 

exponential*sinusoid  form 

IPLT  (12)  Plot  option  15-16 

IPLT  -  0  Plot  routine  not  called 

■  1  Plot  of  the  original  data 

and  model  response  obtained 
-  m  same  as  1  except  that  every 
mth  point  of  the  signals  are 
plotted  (e.g.,  with  m-2  alternate 
points  are  plotted) 

NNPT  (14)  Number  of  signal  points  read  or  17-20  NPT 

through  simulation.  NNPT  should 
be  greater  than  or  equal  to  NPT 
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YYYY  (FIO.O) 
DT  (F10.0) 
BIAS  (FIO.O) 

ANBIAS  (FIO. 


VAR  (FIO.O) 

CARD  #  5.1-5. 


Blank  (unused) 

21-30 

Sampling  interval 

31-40 

Bias  to  be  added  to  data 

41-50 

This  is  for  use  in  Che  simulation 
modes  (ISIM»1  or  2)  when  it  is 
desired  to  study  the  effect  of 
bias  on  data 

Number  of  points  used  for  a  preli-  51-60 
minary  estimate  of  bias.  That  is, 

NBIAS"  Integer (ANBIAS)  number  of 
points  from  the  right  are  used  to 
find  a  crude  estimate  of  bias;  this 
crude  bias  is  subtracted  from  data. 

This  is  useful  only  when  real  data 
is  analyzed  (ISIM"  -1  or  0) ,  and  ignored 
when  simulation  data  is  generated. 

Variance  of  noise  to  data  61-70 

Must  be  used  only  in  the  simulation 
mode  (ISIM"1  or  2)and  should  be  left 
blank  when  real  data  is  analyzed. 


If  ISIM-  -1  or  0  these  cards  contain 
real  data 

If  ISIM*  1  these  cards  contain  (in 
(5F10.0)  format  the  coefficients  of 
the  z-domain  transfer  function;  first 
denominator  coefficients,  and  on 
succeeding  card(s)  the  numerator 
coefficients. 

If  ISIM"  2  these  cards  contain  (in 
(5F10.0)  format  the  coefficients  of 
the  exponential*sinusold  terms;  one 
such  term  on  each  card.  Each  card 
contains  a)  the  weighting  coefficient 
b)  the  exponent  (real),  c)  the  radian 
frequency  of  the  sinusoid,  and  d)  the 
phase  of  the  sinusoid.  If  the  phase  is 
zero,  the  sinusoid  is  a  sine  wave  with 
zero  value  at  k"0. 

Example:  If  y(t)-  e-t+7e~^t  Cos(2t) 
is  to  be  simulated,  then  ISIM-2,  NC0MP«2 
and  these  cards  are  as  follows 

+1.00000  -1.00000 

+7.00000  -3.00000  +2.00000  +1.57079 


1-50  if 
1-80  if 


Leaving  this 
blank  results 
in  the  use  of 
20%  dta  from 
the  right  for 
crude  bias. 


ISIM— 1 
ISIM-0 
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CARD  0  6 

CARD  #  7 

Variable  Name 
(Format) 

IPR  (12) 


Subtitle  card.  This  alphanumeric  subtitle  Is  reproduced 
in  the  output 

Second  option  card.  This  numeric  information  is  also 
reproduced  on  the  output  (for  user  convenience) 


Columns 


IZTS  (12) 


IREM  (12) 


Description 

Print  option.  Increasing  value 
results  in  more  printing.  Use 
0,  1  or  2  for  normal  use. 


Values  3,  4  and  5  are  useful 
for  diagnostic  purposes,  or  when 
resut-s  of  intermediate  computations 
are  needed  (for  example,  the  first-order 
filtered  signals,  i.e.,  the  information 
signals  are  printed  only  if  IPR.GE.4) 

z-domain  to  s-domain  conversion  3-4 

option. 

IZTS  *  -1  conversion  to  s-domain  is 
not  performed 

-  0  conversion  performed;  only 

the  poles  in  s-domain  printed 
*  1  conversion  performed;  in 

addition  to  poles  in  s-domain 
the  s-domain  denominator  and 
numerator  also  printed. 

Number  of  coefficients  in  the  5-6 

numerator  (of  the  z-domain  model) , 
counted  from  the  right,  and  set  to 
zero 

Method  option  5-6 

ISPN  *  -I  pencil-of-functions  method 
employed 

■  1  pencil-of-functions  method 

employed;  noise  added  to 
simulated  data.  Must  be 
used  only  when  ISIM*1  or  2 

■  0  Analysis  of  noise  only 
>-2  Covaria  ce  equations  used 

(for  LPC  or  Prony) 

■-3  Autocorrelation  equations 
used  (for  LPC  or  Prony) 


Preferred/Default 

Value 


Warning:  Do  not  use 


For  LPC  set  IREM- I 9 


For  LPC  set  IREM- 19 
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IF IX  (12) 


9-10 


NFIX  (12) 

I BIAS  (12) 


IB0  (12) 


Noise  correction  option 
(under  development  -  ignore) 

Auxiliary  parameter  for  use 
with  IFIX  -  ignore 

Bias  extraction  option 

IBIAS  =  0  Bias  extraction  not 
exercised 

-  1  bias  extracted 


11-12 

13-14 


Warning:  Do  not  use 
values  other  than  0  or 


Direct  transmission  term  option  15-16 

IB0  ■  0  Constrains  b..*0  in 

numerator  determination 
IB0  -  1  Model  assumes  direct 

transmission  is  present, 
and  is  determined  together 
with  other  coefficients 


MNPT  (14)  Controls  the  number  of  points  used  17-20 

in  the  computation  of  error,  and 
for  printing  (when  IPR.GE.2)  and 
plotting 

QI  (F10.0)  z-domain  pole  of  the  first  order  21-30 

filter  (i.e.,'the  parameter  q  of 
Section  I) 


NPT 


suggest 

QI-0.8 


DFAC  (F10.0)  Auxiliary  variable  for  use  with 
IFIX  (under  development) 


31-40 


SECTION  V 


APPLICATION  EXAMPLES 


Two  examples  will  be  presented  In  this  section.  The  first  deals 
with  a  simulated  noisy  signal,  x(k)«  y(k)  +  w(k)  where  y(k)  is  the  response 
of  a  third  order  transfer  function  and  w(k)  is  a  zero  mean  noise  component. 

The  performance  of  the  pencil-of-functions  method,  with  and  without  Gram 
matrix  enhancement,  is  compared  with  that  of  other  methods  [l]-[3].  The  second 
example  pertains  to  the  transient  response  of  a  conducting  pipe  tested 
at  the  ATHAMAS-I  EMP  simulator.  These  examples  demonstrate  the  effective¬ 
ness  of  the  pencil-of-f unctions  method  as  a  practical  modeling  technique. 
Example  1 
Let 


y(k) 


1  -  1.92z~l  +  z~2 _ 

1  -  2.68a”1  +  2.476z“2  -  0.782z"3 


where  the  sequence  y  is  truncated  at  K«99.  The  signal  under  test  is 
obtained  as 


x(k)  -  y(k)  +  w(k) 

where  w(k)  is  a  zero  mean,  uncorrelated  sequence  with  standard  deviation 
equal  to  0.0316.  The  energies  of  the  signal  and  noise  components  are  1.82 
and  0. 1, respectively ,  hence  the  sequence  under  test  has  a  slgnal-to-noise 
ratio  of  12.6  dB.  The  numerical  data  and  a  plot  of  the  signal  under  test 
are  given  in  Appendix  B. 

The  signal  under  test  was  analyzed  by  the  following  methods  [2]. 

1.  Pencil-of-functions  method 

a)  without  enhancement  of  the  Gram  matrix 

b)  with  enhancement  of  the  Gram  matrix 

2.  Linear  predicive  coder  (LPC) 

a)  using  autocorrelation  equations 

b)  using  covariance  equations 
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3.  Prony's  method  [11] 

a)  using  autocorrelation  equations 

b)  using  covariance  equations  [3] 

The  computer  output  for  these  six  runs  is  given  in  Appendix  B.  Here 

we  summarize  the  fractional  energy  error,  defined  as 

K  2 

2  <yOc)  -  yOO) 

k*0 

V  -  —  -—K -  (28) 

2  y2(k) 

k-0 

A 

where  y(k)  is  the  impulse  response  of  the  model.  Note  that  in  simulation 
mode  the  model  response  can  be,  and  is,  compared  with  the  true  signal. 
When  real  data  is  tested  the  true  signal  is  not  available,  hence  x(k) 
must  be  used  in  this  equation  instead  of  y(k). 

For  the  present  example  the  true  signal  is  available,  hence  the 
fractional  energy  error  is  computed  by  exact  application  of  (28).  The 
values  for  the  various  methods  are  tabulated  in  Table  1. 

TABLE  1 

_ Method _ Fractional  energy  error 


Pencil-of-fuiictions  0.0167 

Pencil-of-f unctions  with  enhancement  0.0016 

Linear  Predictor  (autocorrelation  eqn)  0.1444 

Linear  Predictor  (covariance  eqn)  0. 1567 

Prony  (autocorrelation  eqn)  0. 1431 

Prony  (covariance  eqn)  0. 1397 


The  model  responses  are  compared  with  the  true  signal  y(k)  in 
Figures  6  through  11.  The  reader  is  cautioned  that  the  solid  line  in 
each  of  these  figures  represents  the  true  signal,  and  not  the  noisy 
signal  under  test.  The  latter  is  shown  in  Fig.  B2  in  Appendix  B. 


The  dotted  line  in  each  of  these  figures  represents  the  model  response. 
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true  Impulse  response 
model  Impulse  response 
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Flo.  7.  Comoerlson  of  true  and  model  Impulse  responses 
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true  impulse  response 
model  Impulse  response 


This  example  demonstrates  the  superiority  of  the  pencil-of-functions 
method  over  two  other  widely  used  methods  (for  modeling  impulse  responses) , 
namely  the  all  pole  linear  predictor  and  the  prony  method. 

Example  2 

As  a  real  world  application  we  consider  the  use  of  pencil-of-functions 
method  to  the  transient  response  of  a  conducting  pipe  tested  at  the  ATHAMAS-I 
DIP  simulator.  The  conducting  pipe  is  10m  long  and  1  m  in  diameter.  Hence, 
the  true  resonance  of  the  pipe  is  expected  to  be  in  the  neighborhood  of  14  MHz 
MHz.  Also,  the  pipe  has  been  excited  in  such  a  way  that  it  is  reasonable 
to  expect  only  odd  harmonics  at  the  scattered  fields.  The  data  measured 
are  the  Integral  of  the  E-field;  t.e.,  the  measured  variable  is  a  voltage. 

The  transient  response  used  for  analysis  is  shown  in  Fig.  12  by  the  solid 
line.  The  results  of  analysis  by  the  pencil-of-functions  method  are  given 
in  Appendix  C  for  the  case  of  an  8th  order  model;  the  model  response,  with  an 
error  of  0.0125,  is  shown  in  Fig.  12  by  the  dotted  line.  Model  poles  are: 


fundamental 


-4.280  +  J  67.686  Mrad/s  (-10.794  MHz) 


3rd  harmonic  -22.470  +  J218.200 
curve-fit  pair  -2.543  +  j  12.890 
curve-fit  pair  -16.547  +  j  88.981 


(-34.911  MHz) 
(-  2.051  MHz) 
(-14.404  MHz) 


Mote  that  a  pole  at  the  origin  (due  to  the  integrator)  has  not 
been  obtained  because  we  have  used  the  bias  extraction  option.  On  the 
other  hand,  the  curve-fit  pole  pair  arises  because  the  data  do  not  truly 
pertain  to  the  impulse  response  of  a  finite  order  (lumped)  linear  system. 

Next  the  data  were  differentiated  (actually  differenced),  and 
analyzed  by  the  pencil-of-functions  method.  The  results  are  given  in 


Appendix  C  for  8th  order  analysis;  the  model  response,  with  a  fractional 
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energy  error  of  0.0369  is  shown  in  Fig  13  by  the  dotted  line  (the  solid 
line  depicts  the  differentiated  data).  The  model  poles  are: 

fundamental  -9.453  +  j  71.596  Mrad/s  (-11.494  MHz) 

3rd  harmonic  -25.996  +  j 222. 340  "  (-35.628  MHz) 

5th  harmonic  -113.303  +  J617. 095  "  (-99.855  MHz) 

curve-fit  pair  -22.268  +  j  59.213  "  (-10.068  MHz) 

Here  again  a  curve-fit  pole  pair  arises  because  the  data  do  not 
truly  pertain  to  the  impulse  response  of  a  finite  order  (lumped)  linear 
system. 

We  also  give  below  the  poles  arising  from  a  6th  order  model 
(the  computer  output  is  not  given,  nor  the  graphical  display  of  the  model 
response) : 

fundamental  -9.575  +  j  75. 106  Mrad/s  (-12.050  MHz) 

3rd  harmonic  -14.866  +  J221. 693  "  (-35.363  MHz) 

curve-fit  pair  -14.123  +  J  34.045  "  (-  5.866  MHz) 

The  fractional  energy  error  in  modeling  is  0.0566. 

Here  again  a  curve-fit  pole  pair  arises  because  the  data  do  not 
truly  pertain  to  the  impulse  response  of  a  finite  order  (lump  ed)  linear 


system. 


CONDUCTING  PIPE  TEST 


Fig.  12.  Comparison  of  measured  data  (of  the  response  of 
of  a  conducting  pipe)  and  model  Impulse  response 
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CONDUCTING  PIPE  TEST 


Derivative  of 
Measured  data 

Model  impulse  response 


Method:  Pencil -of-functions 
Enhancement  used 
Bias  extracted 


l 
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Comparison  of  the  derivative  of  measured  data 
and  the  corresponding  model  impulse  response 
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APPENDIX  A 


DETAILS  OF  PROGRAM  POF-FILTER 

LSITINGS 

PURPOSE,  EQUATIONS 
FLOWCHART 
VARIABLES 

FOR  THE  MAIN  AND  SUBROUTINES 
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MAIN  PROGRAM 


Initialization 
Read  option  cards 


Read  or  Generate  signal  y 


listing  of  POF- FILTER 


PROGRAM  "PGF-FILTER" 

IMPULSE-RESPONSE  MODELING 
BY  PtNCIL-OF-FLNCTIONS  M£ThOO 
APRIl  15:0 

*  DtCOUPLtO  QENOM.  ANO  NLM.  DETERMINATION* 

*  FIRST-CROER  FILTERS  USED* 

*  NOISE  CORRECTION  OPTION* 

*  BIAS  CORRECTION  OPTION* 

*  OIRECT  TRANSMISSION  OPTION* 

*  RESULTS  IN  BOTH  Z-  AND  S-  DOMAINS* 

“POF-FILTE  R‘*  MODELS  IMPULSE  RESPONSE  CF  A  SC  ATT  E  RE  R /CHANNEL /NET  WORK 
IT  CAN  Bt  L'StO  IN  SIMULATION  MOOE 
OR  ON  EXPERIMENTALLY  RECOROEO  RESPONSES. 

FOR  COMPARISON  PURPOSES  IT  ALSO  PROVIDES 
ON  AN  OPTIONAL  BA  SIS  THE  FOLLOWING  METHODS 

*  LINEAR  PREDICTOR  (COV  OR  ALTOCORR) 

«  PRONY 

*  LEAST-SQUARES  PRONY 


NOTE!  1350  LINES  OF  CODE.  THIS  CAN  St  REDUCED  SUSTANTIAlLY 
FOR  PARTICULAR  APPl  IC  A  TI ONS  .  IN  PARTICULAR,  ROUTINES  "PLOP"  U7  )  , 
" ZTGS" <157 )  ,  ANO  "FClRT"<EQ11  MAY  3t  ELIMINATED 
IF  ONlY  z-ocmain  results  neeoed 

DIMENSION  F  (2  50  0)  ,U< 80  0) , LU  (8  00 )  ,X  ( 8  00, 11 ) ,  G <  11,  ll')  ,AM(ll,il> 
DIMENSION  GN(11,11),GEST(11,11)  ,GOJM(il,ll)  ,t  111,  1 1 ) ,  £N  (1 1. 11) 
DIMENSION  V<2  2)  ,  VV  (22)  ,A MF  ( 11 1  , SR  <  11 )  ,S I <  11 )  •  SPH  <  1 1 1 
DIMENSION  TITLE(70),iauF(512) ,IDUM(10> 

DOUBLE  PRECISION  OT, AC, 30, ERROR 

COMMON  /OAO/ISPN,OELTA,SIG2,OT,QI,3IAS,  19  IAS,QFAC 

COMMON  /DAI /F  oAP.,t3AR,FESLM,EESLM 

COMMON  /IO/IR,IlT,IPR, IT  PR , IZ PR , IROUNO , IP LT 

RE  H IN05 

MAXPL=800 

M  A  X  3 11 

MAX2=2*MAX 

IR*5 

ILT  =  6 

I S  K IP®  0 

CALl  VE3UAT(MAXPl,U,F,0,10) 

CAll  VcQLAT  (MAXPL,L'U,F,0 ,11> 

CAll  VEQLAT (MAX2, V,VV,0,0> 

WRITE (ILT,2) 

REA0(IR,8) (TITLE (I>, 1*1* 70) 

WR  ITS  (  IlT,  1  c)  (TITLE  (I),  I  *1,70) 

RE  AO (I fi ,8)  (TITLE(I) ,1*1,70) 

RE  AO  ( I C, a) (TITLE (I) ,1*1,70) 

R£AO(IR,4)NPT ,IXX,IYY,N,ISIM, NCCMP ,IPlT,NNPT, 

♦  YYYY,DT, 21  AS,  AN 31  AS,  VAR 
NP1=N*1 
NP  2  *NP  1* 1 
NP3®N*  3 
NPNP2*N+N+2 
npnPi*n^n»i 
IF  (NNPT.EQ.  C)  NNPTSNPT 
Ir (OT.EQ.J .0) DT*1.0 
IFuN*2 


37  (mn-2) 


IF  ( IPlT.  LT  •  u)  IFUN*1 
IF  ( IPuT •  LT  •  G)  IPLTa-IPLT 
FMAX=1.Q 

IF ( ISIM.SQ. 3) GO  TO  63 
IF(ISIM.NE.-l)GO  TO  56 
IMAX=0 

00  58  1=1,  NNPT  » 10 

READ  <IR,  161  XIOUM  <J>  ,J=1,10) 

00  53  J=l,10 

IF(IA35(IDUM<  J)  )  •  GT.  IMA X)  IMAX=IA8S  <  IOUM  (J  )) 

<=I-1+J 
F ( K ) =1 DUM ( J ) 

CONTINUE 
FMAXsIMAX 
GO  TO  45 
CONTINUE 

IF  ( ISIM, NE . 0  GO  TO  49 
R£AD(IR,163)(F(K)  «K=1«  NNPT) 

FMAXxO .0 
DO  41  K  =  1,NPT 

IF  (A6S  (F  <K  I  ).  GT.FMAXIFMaX=A8S  IF  IK} } 

CONTINUE 
NS  I AS=AN3 IA  S 

IF (NSIAS.EQ.O )N8IAS=0 .2*NNPT 

N0I£=NN?T*i-N8IAS 

F3=0.0 

00  57  K=NOlE, NNPT 
F3=F3*F(K) 

F3=F8/N3IAS 
DO  56  Ksl.NNPT 
F(K)s(F(K)-F8)/FMAX 
GO  TO  61 
CONTINUE 

IF  (ISIM.EQ.l)  RE  AO  (IRt  160  )  (V  <1 1  »  1*1*  NP1) 

IF  (ISIM.EQ.l)  READ  <IR, 160)  (V  (I )  ,I=NP2,NPNP  2) 

CAL*.  VcQUAT  (NP1  ,V  (NP2I  ,-1.0,0  ,3) 

IF  <ISIM.EQ.il  CALL  R£SPON  <  F  ,  U  ,  N,  V,  VV  t  NNPT) 

IF (ISIM.EQ.l) GO  TO  61 
DO  60  I=l,NCOMP 

RE  AO  <IR,5>  AMP  (I), SR (I ) ,SI (I), SPH<I) 

WRITE  (IuT,  11)  I, AMP  ( I )  t  S  R  ( I )  ,SI(I)»SPH(I) 

C AUu  SIGNAL  (F ,NNPT,AM?,SR,SI,SPM,OT ,NC0MP) 

CONTINUE 

IFdPuT.Gc.  1)  CALL  PLOTS (IBUF,  512,9) 

REAC  (IR, 8) (TITLE (I) ,1=1,70) 

IF (EOF (IR) • NE • 0 ) G 0  TO  998 
WRI T£ ( IuT • 3 ) 

WR ITE  (  IlT,  16)  (TIT LE(I)  ,1  =  1,70) 

R£AQ(Ifi,6l (TITL£(I), 1=1,701 
WRITE  (  IlT,  16)  (TITLE(I)  ,1  =  1,70) 

3 ACKSPACE  5 

REA0(IR,6) IPR,IZTS,IR£M,ISPN, IFIX,NFIX, 1= IAS, 130 , MNPT ,Q I, OFAC 
NSTRT=2 

IF ( ISPN.EQ.-2 )N3TRT=NP1 

IF(I$PN.EQ.-3)NSTRT-I 

IZPR=0 

IF ( IPR.GE.2 J. ANC. IPR.LE.36) IZPR=(IPR-10 1/10 
ITPSsg 

IF  (IPR.GE.90), ITPR  =  1 
IPR  =  I  =  R-1C-*  (IPR/10) 


38  (m::i~z) 


IF  (ISPN.ll.-2 .AMO.IR£M.iT.NIIS£M=N 

IF (DFAC.EU. 01 DFAC*Q. 1 

133-1 

IF (I3IAS.EQ.G )I3S*0 

IF(I£IAS.lc.0)IBIAS=0 

IROLNO*Q 


CORRUPT  SIGNAl  IF  OESIREO. 

PROCESS  WITH  FI  AST  ORDER  FILTERS 

CONTINUE 

IF (WAR. GE. 1.0 E- 6) CALL  CORLPT (F , X, VA R , NPT , HAXPL) 
CONTINUE 

IF (VAR.GE.l.QE-6) GO  TO  95 
00  30  K  =  1,NPT 
X  ( K  »  1)  =F«)  +5  IAS 
CONTINUE 
I NT--1 

IF ( ISPN. L£ • -2 ) I  NT  =2 

IF (NP1.GT. 1)C ALL  FILTER ( X  «NFT  ,NPi, MAXPu , I NT l 

COMPUTE  GRAM  MATRIX 

NPP*N?1 

IF ( I3IAS.NE.0 )NPP=NP2 
DO  44  1*1, NPP 
DO  44  J=1,NPP 
A 0 - Q  .u 

IF (ISPN. EQ.O.ANO. IROUNO.EQ. 0)  GOTO  43 
DO  42  K=NSTRT  » NPT 
AO=AO*X (K» I I* X (K, J ) 

GN(I,J»=AO*DT 

IJ1-IABS(I-J)+1 

IF(ISPN.EQ.-3.AND.I.GE.2)GN(I,J)-GN(1,IJ1) 
GOUM(I,J)*GMI,  J) 

continue 

IF (IROUNO.EQ. 0) G(I,J)*GN ( I, Jl 
CONTINUE 

IF(ISPN.NE.O.OR.IROUNO.NE.O) 

1C  ALl.  GKROCT(GN,E, CET , V,NFP,NPP, MAX, 1) 

IF (IROUNO.EQ.  0) WRITE (IlT,171) OET 
IF ( IROUNO.EQ. 1) WRITE (ILT.172) OET 
IF (IPR.GE.l)CAtL  PRTMAT(GN,NPF,NPP,MA  X, -1 ) 

WRITE ( IwT ,  1) 

I RQ*IRCUNQ 

IF  (IROUNO.EQ.  OHROUNO=IROLNO>1 
IF (IRO.EG.O.AND.ISPN.GT.-l)GO  TO  410 
IF (IFIX, EQ. -1 )  GO  TO  203 


ESTIMATE  OF  ♦*  G 

CALL  3UIL0Z(AM, V,NPl,NPT ,MAX,NFIX) 

- - NF1  REPLACED  3Y  NPP  NEXT  3  CARDS— 

CAll  FIX(G0UM,AM,G£ST,E,V,NFP,NPP,S  IG2.MA  X,  IFIX) 
IF(IFIX.EQ.I) WRITE (ILT ,462) SIG2 
CAL-  GKROCT  (GES  T,  E  ,OET ,  V , NPP,  NPP, MA  X ,  1) 
WRITE(IlT,162)DET 

IF ( IPR. GE. 1)C ALL  PRTM«T(GEST,NPl,N®l,MAX, 0) 

DC  15-,  1*1,  NP1 
DO  154  J  *  1 «  NP  1 
GDUM  (I  ,  J)sGtST (  I,  J) 
nfix*nfix-i 


29  (MZ'J-4) 


IF  (NFlX.Gc.il  GO  TO  156 

ISKIP=1 

IR0tN093 


DETERMINE  NUMERATOR 
CONTINUE 

IF ( IBS. EQ. 1)1  BIAS9 1 

IF  ( ISPN.EQ.  G)  GO  TC  596 

CALL  VEQLAT  (NP1,V  (NP2) ,VV,0,Q) 

V (NP2) *-1.0 

CAll  R£SPON(X (1  ,1) ,U,N,V, VV.NPT) 

IF (IPR.G£. 41 WRITE (ILT, 174) 

IF (IPR.Gc.4IW  RITE  ( ILT  ,  210  )  ( X  (  K  ,  1)  ,  K  *1 ,  NPT  > 

CALL  FILTER  (X  ,NPT  , NP1- I fi£M , MA  XPL.2) 

L*N-IREM 

IF(I0IAS.NE.O  )  L  *L*i 

LPlsL+1 

LP2=L*2 

IF  ( IBI AS.NE.O I  CALL  VEQUAT (NPT  .X (It. Pit ,L,  0, 11) 

CALL  VEQLAT  (NPT, X  (1.LP2)  ,F, 0,1) 

CALL  VEQLAT  (NPT  ,X(1,LP2)  .BIAS, 0,4) 

IF(IPR.GE.5)WRIT£  (ILT,  110)  (  (X  (K  ,1)  ,  K*l,  NPT)  ,1*1,  LP2) 

L=L+IS0 

LP1*L*1 

00  216  I*1»L 

DO  216  J*i rLP 1 

AD*Q.O 

DO  215  <=1, NPT 

AO*AO+X(K,in-IBQ  >*X  (K.J+1-I30I 
G(I,J)*AO*OT 

IF(IPR.Gt.5)CALL  PRTMAT (G  ,L  .LPltMAX,  2  05) 

CALL  GKROCT (G  ,£ ,OET, VV ,l ,L ,MA X, 0) 

IF ( IPR.Gt. 5)C ALL  PRTMAT (£ ,L ,L  .MAX, 20  7) 

CALL  VEQLAT  (NPl.VV, AMP, 0,0) 

00  220  1*1, L 
A  0s 0 .0 

DO  219  J  =  1,L 
AO*AOfE(I,J)*G(J,LFl) 

VV  (  I) * AO 
FM£AN=0.0 

IF  (131 AS.NE.O )FM£AN=VV(_ ) 

CALL  VEQUAT  (I  REM,  VV  (N>13  0 -I  fi£  M  +  l)  , A  MP,  0  ,0  ) 

V (HP 2) *0.0 

CALL  VEQLAT (N +130 , V (NP3-IB0 I , VV,0, 1 ) 
IF(IPR.GE.2)MRITE (ILT,12)FMAX 
WRITE (ILT, 303 )F MAX 
WRITE  (Iv.T,  210  )  <V(  I)  ,I*1,NP1> 

WR ITE ( I-T, 210 ) (V(I),I*NP2,NPNP2I 
IF  (IB  I AS.NE.O I  WRITE (IlT, 30  5) FME AN 

model  response,  ano  error. 

IF (MNPT.EQ. 0) MNPT*NNPT 

CAlL  VEQLAT  (N PI  ,V  (NPZ) ,- 1.0, 0,3) 

CAlL  RES  PON  (X  (1  ,2)  ,U,N,V,VV.MNPT) 

CALL  VEQLAT (MNPT, X (1, 2 ) ,FM£ AN , 0 .4) 
c  R ROR* 0  »  u 
FFSlM*0. 0 
00  213  <si , MN FT 

x  (<,d*f  (o*fmax^3Ias  (yAi::-c) 


X  (<,  2)  *X  «,  2)  *FMAX 
FFSLM  =  FFSuM,X  «  ,1)*X  (K  ,1 ) 
X(«,3)*X  « , 1)  -X  (<  ,2) 


13  £RROR  =  £RROR*X  (K,3)*X  (K,3) 

FFSuM*FFSUM*OT 

ERROR=ERROR*OT 

RAT  I0=£RR0R  /FFSUM 

WRITE  (Ii-T, 304 ) ERROR, FFSUM, RAT  10 

IF (IPR.3E,2)WRITE (ILT, 112) 

IF (IPR«GE.2)WRIT£ (ILTtllO) (X(K, 1) ,K *1,MNPT) 
IF  (IPR«GE»2)WRITE  (ILT. 113) 

IF <IPR,GE.2IWRITE  (ILT, 110 )  (X (K , 21 ,K=1,MNPT) 
IF(IFR.GE.2)WRITE (ILT, 114) 


IF (IPR«GE,2)WRIT£ (ILT, 110)  (X (K, 3) ,K*i,MNPT) 

OELT*OT 

IFCIZTS.Gc.  OCALL  ZT05  (V  Cl)  ,V  (NP2)  ,  N,  CELT  ,  I  ZTS) 

IF ( IPW T. EQ , 0) GO  TC  239 
IF(IPlT.EQ.I)  GO  TO  236 
00  230  1*1,2 
<K*0 

DO  230  K»lt  MNPT.IPLT 
KK*t<K*l 

«0  X (KK,I)=X(K,I ) 

mnpt=mnpt/iplt 
>  b  T  0  *  0  •  0 

OELT*OT*IPlT 

call  PLOP ( MNP T, IFUN,X,MAXFL,TO,OELT ,lhY,lHT,IBUF) 

19  CONTINUE 

FORMAT  STATEMENTS 

FORMAT (14,612.14, 6F10.0) 

FORMAT  (5F10. 3) 

FORMAT (812, 14, 4F10.0) 

FORMAT  (70A1) 

FORMAT  (2X.70A1) 

3  FORMAT (1C(5X,  F5.0)  ) 

4  FORMAT (10(5X,  15)) 

FORMAT (2X, 12,  ♦  AMP**fFa,2,* 

1*  PHASE **»Fl  Q,  4) 

FORMAT  (2X,*WA  VEFOFMS  AND  NUMER. 

0  FORMAT  (10X,bhG  MATRIX) 

0  FORMAT  (10X.6HM  MATRIX) 

2  FORMAT  (1QX,  11HGEST  MATRIX,*  (OET=* ,G13 . 6 »♦ >• > 

0  FORMAT  (2X,5  (2X,G13.6)) 

10  FORMAT (20(1X,F5. 2)) 

0  FORMAT (2X,10(F7. 4) ) 

2  FORMAT  (2X,*0RIGINAL  SIGNAL  ( INCLUDES  BIAS,  IF  ANY)*,/) 

3  FORMAT (2X,*IMPL.RESP  OF  MCOEL (INC. 3 -HAT ,1 F  I6IAS.NE .0)*,/ ) 

4  FORMAT  (2X,*cRROR*F«)-FR£C(K)  *,/) 

6  FORMAT  (10X,  14HNOISY  X  MATRIX) 

5  FORMAT  (10X, 6HX  MATRIX) 

1  FORMAT  (10  X ,  16  mTRUE  GRAM  MATRIX,*  (0£T=*  ,G  13 .6,*)  *) 

2  FORMAT (10X, 17hNOlSY  GRAM  MATRIX,*  ( OET»*, G1 3. 6,* )*) 

4  FORMAT  (2X,*IMPL'LSE  RESPONSE  OF  1/A(Z)*) 

3  FORMAT C5F1Q.0  ) 
i  F  CRMAT (1015) 

3  FORMAT  (10F3. 6) 

3  FORMAT  (2X,  *tST 
♦  FORMAT (/,2X,*SS 
**RATI0**,G13.6,/> 

5  FORMAT (2X, *ES TI MATEO  M£A N** ,G 13 . 5) 

?  *CRMAT(*  ESTIMATED  NOISE  VAR**,G12.5) 

FORMAT  (/)  41  (yji::-6) 


S**,Fl  0,4,*  J*,F  10,4, 


SCALEO  3 Y  XMAX=*,F12.  5) 


TF  3  (Z)  /A  ( Z  *  *  ,*  (FMAX**,E12.5,*)J 

ERR0R**»G13»6,*SS  S IGNAL* *, G13.6, 


) 


FORMAT (1H1) 

FORMAT (////) 

GO  TO  ill! 

CONTINUE 

CAu.  PLOT  ( 0  .*  0.  *999) 

STOP 

ENO 
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SUBROUTINE: 


BUILD R 


PURPOSE:  To  generate  a  conversion  matrix  to  go  from  A) 

to  parameters  a  of  A(z).  See  equation  (9)  1 

of  Section  I 

EQUATIONS: 

Conversion  matrix  (shown  below  for  n«3) 

~1 

-3  1 

3  -2  1 

-1  1-1  1 


3  2 

-  DiagCq  ,  q  ,  q,  1} 


ROUTINE  VARIABLES 

A  Conversion  matrix  (not  to  be  confused 
with  the  denominator  polynomial) 

X  Vector  which  brings  in  [ A>^  ....  /Dn 

to  the  routine  and  takes  back 

the  denominator  parameters  [aQ  ....  a^ 

N  System  order  plus  one 
MAX  Maximum  dimensionality  of  matrix  A 


FURTHER  DESCRIPTION: 
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SUBROUTINE  3U ILOR  (A,X«N»MAX) 


CONVERSION  MATRIX*  REVERSE  FOF  PROCESSING  —  I.R.  MODELING 
DIMENSION  A  (MAX  *1)  ,  X  (l)tY  (20) 

DOUBLE  PRECISION  OT,Y 

COMMON  /DAO /I SPN,OELTA,SIG2,OT,QI,3 IAS, 15 IAS,DFAC 
COMMON  /I0/IR,ILT,IPR, IT  PR , IZPR , IR3UND , IPLT 
N  Ml *N- 1 
00  11  1*1, N 
Y (11=0.0 
DO  11  J=1,N 
A (I, J) *0.0 
A (N, N) *1.0 


DO  20  J J*1 , NM  1 
J*N-JJ 

DO  15  KK  =1 , 2 
K*<K-1 

DO  15  I=J,NM1 

A  ( I +K, J) =A  (I*K,J»  ♦AII+1,  J*l)*(1.0-K-K> 
CONTINUE 
00=1.0 

CHANGED  Thru  22  3-17-30 

DO  22  1 1  * 1 , NM  1 

I  *  N  - 1 1 

QQ*CQ*CI 

DO  22  J*l, N 

A (I, J) =QQ* A  (I  .J) 


DO  25  1*1, N 

IF(IPR.GE.3)HRITE (ILT.5) (A (I, JJ) ,J J =1 ,N> 
00  25  J=1,N 
Y  (I)*Y  (I )  ♦  A  (I  ,J)*X(J) 

DO  26  1=1, N 
xm=Y(I)/YCl> 


IF(IPR.GE.3)wRITE  (6,7)  (X(I)  ,1*1, N) 

FORMAT  (2X,  10G12.5-) 

FORMAT  (*  ESTIMATED  PARAMETER  VECTO  R*  ,  /  ,  1  CG  13 . 6) 

RETURN 

END 


SUBROUTINE: 

PURPOSE: 


EQUATIONS: 


BUILDZ 

To  calculate  unit  noise  covariance  matrix 
for  reverse-time  first-order  filtering  case 
(under  further  development) 


sv  — 

Z  -  I  (K-k+1)  m(k)  mT(k) 
k-1 

where  m(k)  is  the  vector  of  unit-pulse  responses 
of  the  connection  filters 


ROUTINE  VARIABLES 

Z 

R 

NP1 

NDIM 

NPT 

NFIX 


Covariance  matrix  for  unit  noise 

Work  vector 

system  order  plus  one 

Maximum  dimension  of  the  matrix  Z 

Number  of  points  in  signal 

Not  used  (blank) 


FURTHER  DESCRIPTION: 


siasOuTi.Nc  at  iuoz  <z,r  .n?i, nft  ,nuim, nf  ix> 


AO  TE  SnAT IVE  NOISE  CCV  PGM  FOP  *GNN * 

DIMENSION  Z  (NOIM,  1 )  ,R  ( 1) 

D0U30E  PRECISION  OT 

COMMON  /OAO/IS?N,OtOTA,SIG2,OT,CI,3IAS,  13  IA  S,OFAC 
COMMON  /IO/IR  ,IuT, IPR, ITPS, IZFRfIROUND, IPLT 
Q  =  QI 

T I H£*DT*NPT 
IOPT=NFIX*l 

GO  TO<201, 101,101,201) , ICPT 

SC*OT 

NSNP1-1 

R (1) =1.0 

DO  1  1*1, NPl 

IF  ( I .Gt • 2) R  (I  )*R(I-1) 

DO  1  J=l, NPl 
Z (I,J) =0.0 
DO  2  <=1,NPT 
NPTK=NPr*l-< 

DO  3  J*1 *N?1 
DO  3  I *J,NP1 

Z(I,J)=Z(I,J)  ♦R(I)  *R(J)*NPT< 

P ( 1 ) *0 . 0 
DO  4  I = 1 , N 
RtMUWlIUl  ♦R(I) 

CONTINUE 
DO  7  J=1,NP1 
DO  7  I = J , NPl 

Z  <1,  J)=Z(I.  J) *DT** (I ♦ J-l 1 

GO  TO  301 

CONTINUE 

DO  210  J  =  1 , nP  1 

OO  210  I =  J , NP 1 

IF(I.EQ.1)Z(1,1)=TIME 

Z  <I,J)*<TIM£*MI+J-1>  )/(l4j-l) 

HR  I  T£  (  IlT,  1 61 ) 

CONTINUE 
DO  174  1*1, NPl 
OO  163  J  =  I  ,  NP  1 
Z<I,J)*Z(J,I) 

IF(IPR.G£.2)WRITE (IUT«22G) CZ(I»J),J=i,NPl> 
CONTINUE 

FORMAT  <10X , •  QUANT.  NOISE  *) 

FORMAT (10X, *3  IAS  EFFECT*) 

FORMAT  C2X,5  (2X.G13.6)) 

RETURN 

END 
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SUBROUTINE: 


CORUPT 


PURPOSE:  This  routine,  useful  in  simulation  mode,  can  be 

called  when  it  is  desired  to  add  random  noise 
of  given  variance  to  the  signal 


EQUATIONS: 


ROUTINE  VARIABLES 

F  Input  signal 

X  First  column  of  X  would  contain  the  corrupted 
signal;  the  second  column  temporarily 
contains  the  noise  added  to  signal 

VAR  Variance  of  noise  added  to' signal 

NPT  Number  of  signal  points 


NDIM  Maximum  column  dimensionality  of  X 


FURTHER  DESCRIPTION: 


This  routine  needs  a  library  routine  to  produce 
random  numbers  (gaussian,  zero  mean  and  uncorrelated) 


SU3  RQU  TINE  COftUPT  (F,X,VAR,NFT ,NOlM) 


ADOS  NO IS; 

DIMENSION  F  (1 )  , X (NQIM, 1) 

OOU3l£  PRECISION  DT,AQ,30 

COMMON  /DAO /I SPN,0£LTm ,SIG2,0T,QI, 3  IAS, 13 IAS,DFAC 
COMMON  /OA1/F  £A  R ,  E  3AR  ,  FE  SU  M  ,E  E  SUM 
COMMON  /IO/IR,lLT, IPR, ITFfi,IZPR,IR3UN0,IPLT 
F  3  A  R  =  0  . 

E 3  ARsQ • 

FESUM=0. 

£ESUM=0. 

WRITE(lLT»48S)VAR 
IS=2458169 
IS2=397665 
SIGMA  =  SQRT  (  VA  R ) 

CALL  NRML(NPT  ,1,1,0. , SIGMA, IS,  IS2,X (1,2), 0) 

00  26  K=1,NPT 
X  (K,l)sF  (K)  +3IAS>X(K,2) 

00  211  K=1,NPT 
F3  =F (< ) ♦6IAS 
FSAR=FeAR*F3 
E9AR=£3AR+X  (K  ,2) 

££SLM=£ESUM+X  (<  ,2)  ♦X  (K,2  ) 

FESUM=FESuM+F S*X(K,2) 

IF  (ISPN.EQ.  0)  X(K,  1>  =  X(K,2) 


CONTINUE 

EESUM-E£SUM*OT 

FESUM=2.0*FESUM*OT 

FoAR=F3AR/NPT 

EBAR=£SAR/NPT 

WRITE (ILT, 4c2)FeAR,EBAR,FESUM,E£SUM 
IF (IPR.Lt. 2)  GO  TC  411 
WRITE ( ILT , 8 ) 

WRITE  (IlT, 110  )  (X(K,1) ,<=1,N FT ) 

IF(ISPN.EQ.O) GO  TO  411 
WRITE  (ILT, 18) 

WRITE  (IlT,  115  HX  (K,2)  ,K*1,NPT> 

WRITE(IlT,1) 

CONTINUE 

CONTINUE 

FORMAT  STATEMENTS 

FORMAT  (10X,  16HR0UN0E0  F  SIGNAL) 

FORMAT  (1QX, 16HR0UND0FF  cRROR  E) 

FORMAT  (2X,  5  (2X,G11.4)) 

FORMAT (20 (IX, F5.2) ) 

FORMAT  (IX, 20(1X,F5.3) ) 

FORMAT  (2X,5hF5AR=,  Ell.  4, 6h  E9  AR  =  ,£  1 1 . 4,  5H  F  £2  =  ,E  11 . 4,  4h  EE 
FORMAT  (2X,*VARIANC£  OF  N  0  ISE  =  * ,  E12 . 4) 

FORMAT  (/> 

RETURN 

END 
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SUBROUTINE:  filter 

PURPOSE:  To  produce  the  information  signals.  Specifically, 

this  routine  performs  reverse-time  first  order 
filtering  upon  the  given  signal  (stored  in  the 
first  column  of  the  matrix  X) 

EQUATIONS: 

X(k, i+1)  -  q  X(k+l,i+l)  +  X(k,i) 


ROUTINE  VARIABLES 

X  Matrix  of  information  signals.  First  column 
brings  in  the  signal  to  be  processed 

NPT  Number  of  signal  points 

NP1  Number  of  information  signals  (model  order  plus  one) 

NDIM  Maximum  column  dimensionality  of  X 

INT  Option  parameter.  -1  for  reverse-time 

filtering,  +2  for  unit-shifts  (or  delays) 


FURTHER  DESCRIPTION: 
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LiRC^TlN-  F I  u  T  E  R  (  X  ,  N  P  T  »  N  P 1  *  N  D I M  ,  I N  T  ) 


DIMENSION  X(NDIM,1) 

DOUBLE  PP.E C  IS  10  N  0T,SC,30 

COMMON  /DAQ/ISPN,D£UTA,SIG2,0T,GI, 3  IAS, 12  IAS, OFAC 
COMMON  /IO/IR,ILT,IPR,  IT  PR , IZ =R , I ROUND, IP LT 

GENERATE  FIRST-ORDER  FITE  R  FRCCESSEO  SIGNAlS  FROM  n*TA  IN  X(K,1 > 
I N  T  =  1  C»  3  FOR  FORWARD,  -1  FOR  REVERSE  TIME 
FIRST-ORDER  FILTERING 
I N  T  s  2  FOR  UNIT  CElAYS  (X  <K  ,  If  1>  =X  << -1  ,  I  >) 

N*NPi-l 
NP2=NP1+ 1 
IOPTslNT  f 2 

GO  T0<51,11,11,91)  ,IOPT 

FORWARD  FIRST -OROE R  FILTERING 

CONTINLE 

DO  40  J  =  1 »  N 

JJ=Jfl 

X  (1,JJ»=X(1,1) 

DO  40  K=2,NPT 
Kl=<-1 

X (<,JJ)=GI*X(K1,JJ)*X(K,J> 

CONTINUE 
GO  TO  70 

REVERSE-TIME  FIRST-ORDER  FILTERING 

CONTINUE 

DO  60  J  =  l«  N 

jjsjfl 

X  (NPT,  JJ)  =X  (N  FT  ,1> 

X (NPT , JJ ) =0 .0 
3C*X(NPTtJJ» 

DO  60  KK=2,NPT 
<snFT+1-K< 

<- 1  s  K  ♦  1 

CHANGED  NEXT  CARO  3/17/80 
30=30fQI*X(K, j) 

30  =  GI*3QfX  <K,  J) 


X  <K, JJ) =EO 
CONTINUE 

IFdSIAS.EQ.O  IGO  TO  62 
IPWR  =  I3IAS-  1 
DO  61  K  '<  =  1 »  NP  T 
T IME=OT*KK 

<  =  N  FTf 1- << 

X  (K,NP2)  =T  I  ME  *•*  IPWR 

CONTINUE 

GO  TO  70 

GENERATE  UNIT  OElAYS 
CONTINUE 
DO  33  1=2, NP1 
11=1-1 
X  <1,I)=3.0 
DO  93  K=2,NPT 
<1=<-1 

X  <<,I) =X  {<1,1 1) 

GO  TO  31 
CONTINUE 
S  C  =  1 . 3 

DO  33  1=2, NP1 
SC=SC*CT 
OC  5  0  <*  1  f 

<  (<  .  I)  =  3C*X  «  ,1 )  SO  (FILTER  1) 


1  CONTINUE 

IF  <IFR.-T.4IGC  TO  99 
IF (IROcNQ.EQ. 1)  WRITE  <IlT, 176)  INT 
IF (IROLNO.EG. G>  WRITE <IIT,179>  INT 
00  190  1  =  1, N?  2 

3C  WRITE  (IUT.liu > <X<K,I) ,K=1,NPT> 

WRITE  ( ILT, 1) 

9  CONTINUE 

112  FORMAT ( 4  ( IX, F12, 6)) 

1C  FORMAT  (10(1X,F7. 3)1 

79  FORMAT <IOX,*F IL TER, INT=* , 12,*  FOF  PROCESSED  X* ) 

7c  FORMAT  <10X,*FIL TER, INT=*,  12,*  NOISY  FOF  PROCESSED  X*) 
FORMAT  (/» 

RETURN 

END 
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SUBROUTINE:  fix 


PURPOSE: 


EQUATIONS : 


To  enhance  the  Gram  matrix  of  the  information  signals 
and  to  effect  noise  corrections 
(Under  further  development) 


where  P  is  the  unit  noise  covariance  matrix 


ROUTINE  VARIABLES 

G  Noisy  Gram  matrix  of  reverse-time  first-order 
filtered  signals 

P  Covariance  matrix  of  unit  noise  (also  reverse-time 
first-order  filtered) 

C  Corrected  matrix 
D  Work  matrix 


N  Dimensionality  of  the  matrices 
NC  Not  used 

SIG  Estimated  noise  variance 
NDIM  Maximum  dimensionality  of  the  matrices 

IFIX  Option  parameter  (Use  IFIX»1) 


FURTHER  DESCRIPTION: 


Under  further  development 


ESTIMATE  NOISE  INTENSITY  SIG  (ASSUME  WHITE  NOISE) 
CORRECT  NOISY  MATRIX*  C 

IP]  CENOTE3  NOISE  MATRIX  FOR  LN  IT  NOISE 

NC  IS  THE  NCN2ES0  SU5  MATRIX  OF  P  =  C  CV  OF  NOISE 


n  I  MENS  ION  G  (NOIM,  1)  ,P(NOIM,ll  ,C  (NOIM,l)  ,0  (NQIM,1)  ,X  (1) 
COMMON  /OAO/ISPN.OELTA  ,  SIG2  ,0  T.QI,  3  IAS,  IE  IAS.QFAC 
COMMCN  /IO/IR,ILT,IPR,ITFR, IZ PR ,1 ROUND, IP LT 
OOUSuE  PRECISION  SUMOET.OT 
SI*SIG 

IF  ( IFIX.EQ. G) GO  TC  51 

G0?=G(1,1)  /P<  1,1) 

jct=o 

S IG=Q •  0 

jct=jct+i 

suMOET*a.a 

CAuL  GKROCT (G  ,0,G0ET,X ,0 ,N,NO IM ,0 ) 

IF  (JCT.£G.l)D£TG=GOET 
00  7  I  * 1 , N 
DO  7  J=1,N 


SUMCET *SL'MOcT  +0  (I  ,J)*P  (I,J) 

IF (StMOET.LT. 0. 0. AND.IFIX.NE.2) GO  TO  11 
SI  =  1. JCJ/SUMOET 

WRITE ( ILT, 32) JCT, ICT , GOE T , SLMOET,S I 
IF  (SI/GOP.GT.  0.1)  WRITE  (ILT,  31) 

ICT=0 
CONTINUE 
00  3  1*1, N 
00  3  J=1,N 

C (I,J) *G  (I ,J) -S  I*P  (I»J) 

IF  (IFIX.cQ. 0) GO  TC  11 

CAul  GKROCT (C ,0 ,COET ,X ,0 , N , NO IM,0) 

IF (COET.uT.O.C.OR.COET.GT.GCET) ICT=ICT+1 

IF  (ICT.GT.5)G0  TO  10 

IF(ICT.GT.0)SI*SI/2.0 

IF  (ICT.GT.O)GO  TO  51 

IF (JCT.GE.5JG0  TO  11 

THR=QcTG*QFAC 

IF  (JCT.EC.DWRITE  (It  T ,  33 )  OETG  ,  OFAC ,  THR 
SIG*SIG*SI 

IF  (COuT.GT.  ThR)  CALL  MEQUAT  (N,N,G,C,N0IM,1 ) 

WRITE ( ILT,34) JCT, ICT,G0£T,SI, COET 
IF  (COET.GT .TriR)  GO  TO  3 

FORMAT  (2X, ♦NOISE  VAR  EXCESSIVE,  SIG/GOP.G  T.  0.  !♦) 
FORMAT (IX, *J, I  GOE  T,  SCMQET,  SI  **  ,212*  4E1 1, 2) 
FORMAT (1X,*GQET,0FAC,ThR:*,6E 11.2) 

FORMAT  (IX,  ♦  J,  I  SI, COET*, 21 2,  4£  11. 2) 

P£  TtfiN 
ENO 
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SUBROUTINE: 


GKRDCT 


PURPOSE:  Basically,  this  routine  finds  the  cofactors  and/or 

the  inverse  of  a  square  matrix.  It  also  calculates 
the  denominator  parameters  through  pencil-of-functions 
method  (if  ISPN.GE.-l)  or  the  LPC/Prony  methods  (ISPN.LE.-2) 

EQUATIONS: 

-f/l  -n 

A(_l)  -  (?*)  Z  v^/ 

c  -  i 


ROUTINE  VARIABLES 

X  Gram  matrix  of  information  signals 

Y  The  adjoint  or  the  inverse  matrix  of  X  is  returned 
in  Y 

DET  The  determinant  of  X  is  returned  in  this  variable 
XLAMDA  Vector  of  computed  denominator  parameters 
NN  Not  used  (blank) 

N  Dimensionality  of  X  and  Y 
MAX  Maximum  dimensionality  of  X  and  Y 

IOPT  Option  parameter;  0  for  finding  the  inverse 
and  determinant  of  the  matrix  X,  1  to  find 
the  denominator  parameters 


FURTHER  DESCRIPTION: 

Scaling  has  been  introduced  in  this  routine  to 
enable  accurate  computations  even  for  high  order 
modeling  (say  6  to  10). 

This  routine  calls  BUILDR  to  go  from  Di  to  the 
parameters  a^^  of  A(z). 
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SU5  SCU  T INE  GKRDCT  (X , Y , DE T , XwA MOA, NN ,N ,MAX , I  OPT) 


DIMENSION  XlAMO  A ( 1) 

dimension  x (max ,n ,r (max.i) 

DOUcLE  PRECISION  A,B,C,D,E 

DIMENSION  NUM  (2,10>,SCA»(10),RSC(10),Z(10,10) 

00U3LE  PRECISION  0T,SC,A0,3D 

COMMON  /OAQ/ISPN,OEUTA,SIG2,OT,QI,3IAS, IBIAS,QFAC 
COMMON  /IO/IR,IwT, IPR, IT F R , IZ PR , IROUNO, IPLT 
IGKRsi 

IF  (N.NE.i) 50  TO  3 

Y  (1,1) *1.0/X(1,1) 

DE  T*X ( 1 , 1) 

GO  TO  61 
CONTINUE 
NMAT=N 
LPC  =  Q 

IF  (ISPN.LE.-2  .AND.IOPT  .£Q.1)LRC=1 
IF  (wFC.EC.DN  f-AT  =  N-l 

- SCAwE— - 

DO  11  I  s  1 ,  N  f*A  T 
SC  Aw  ( I ) =1. 0 

IF  (X(I*L PC, I* L?C) .GT.O. IE-20) SCAL(I) =SQRT  (X (I ♦LPC , I ♦LFC ) ) 

RSC  (I)  sl.O/SCAL  (I ) 

CONTINUE 
DO  6  I *1 ,NMAT 
DO  6  J*1,NMAT 

Y  (J,I)*X  (J*LPC,  I*UPC)*RSC  (I)*RSC(J> 

Z  ( J ,  I )  sY  ( J ,  I) 

IF(ITPR.Gc.l)  CAU-  PRTMAT  (Z,NMAT,NMAT,10,0) 

- BEGIN  GK  REDUCTION - - 

A  =  1 . 0 

00  43  I  =  1,NMAT 
8*0.0 
U  =  I 
M*  I 

FINO  LARGEST  ENTRY  A(U,M)  IN  THE  LOWER  OIAGONAl  SUBMATRIX 

00  Id  J  =  I,NHAT 
00  la  K  =  I  *  N  HA  T 

IF  (A3S  (Y  (<, J) ) .UE.OAcS (3) ) GC  TO  16 
3*A3S( Y (K, J)) 

L  =  < 

Ms  J 

CONTINUE 

Interchange  rows 

IF  (l.EC.DGC  TO  24 
00  23  J  S  i  9  NMA  T 
C  =  Y  (l»  J) 

Y (w,J) *Y (I, J) 

Y (I, J) *C 

INTERCHANGE  CCwUMNS 
1*  (M.-C.DOO  TO  29 

55  ( GKRDCT- 1) 


00  Za  J=1,NMAT 
C=Y  ( J.  M) 

Y ( J.M) =Y  (J,  I) 

6  Y(J,I)=C 

6  NLlMCl,I»=i. 

NUN (2  » I )  =M 
3=Y  (I,  I) 

Y(i,n*A 
00  42  J  =  1,NPAT 
IF(J.£C.I)G0  TO  42 
Ca“Y  (I  *  J) 

YCI,J>*0.0 
00  41  Kal.NMAT 
0*C*Y(Kf I) 

E*3*Y (K.  J) *0 

I F ( 0 ASS (£) . lT .1. 00-10* 0A5S(D) ) £=0 . 0 

1  Y(K,J)a£/A 

2  CONTINUE 

5  A  *  3 

RESTORE  COLUMNS 

DO  56  1*2 .  N  PA  T 

JsnMAT+i-I 

K=NUM(2.J) 

IF(K.£C.J)GO  TO  52 
00  51  L  =  ltNMAT 
C  =  Y  (K*  L) 

Y  <«,L>  *Y  ( J  f L) 

Y(J,U*C 

K  *NUM ( 1  *  J) 

RESTORE  ROWS 

IF(K.cC.J»GO  TO  56 
00  57  L  *  1 . N  MA  T 
C* Y  {_.  <) 

Y  (L.K)*Y(_,J) 

Y  (t i J) *C 
CONTINUE 
OET=A 
DSC Ac=l. 0 

00  55  Ial.NPAT 

OSC A».=  SC Av.  ( I)  «0SCAL 

OETsOET*SCAl( I) *SCAL (I) 

IF  (ITPR.GE.l)  WRITE  (ILT  .3371  OET,A.(RSC  (I) ,  1= 
7  FORMAT (1X»*0E  T.A,  RSC(I)»  *,7E11.4) 

IF CIOPT.EQ.l. ANO. ISPN, GE . - 1 IGO  TO  61 
IF  (ITPR.GE.l)  CAU.  PRTMAT  (  T  .  N,  N  ,  MA  X  ,  0  > 

IF(  ITPR.Gc.  1)  CALL  PSOMAT  (Z,  Y.NMAT,  10,  MAX,  0) 
00  60  I s  1 »  N  PA  T 
00  60  J*lt  N  PA  T 

Y (I  * J) *Y (I .J)  *RSC ( I) *RSC (J) /A 
CONTINUE 

IP(ICPT.NE.l)  GO  TC  1000 
IF (ISPN. LE. -2 )G0  TO  440 
IF ( Y  (1,1) .LT. 0. 0) GO  TO  1000 


SCsi.O 

00  200  I *2  »  NM  AT 
co i cp* pt 

5SC  ( I)  *RSC  < I)  /RSC  (1) 

IF  (IG<S.E3. 0) XLAMCA(I) »RSC(I) "Y (I,l)/Y (1,  1) 
IF  (IG<S.El.  C)  GO  TC  196 

oc  (sys.: 

k 


l.NMAT) 


XT- 2) 


A  =  Y  (i,  i) 

IF<Y  (1,1). LT.  C.  □ )  A  =  0.0 
IF  (IG<R.EQ.2) A=A3S  (Y  (1,1 ) ) 

XLAMOA  (I)=RSC  (I)*OSQRT  tA/Y (1,1) > 

IF(Y(I,1).LT. 0. 0) XLAMOA(I) =-XL AMOA ( I ) 

XLAMOA  (I)=SC*XLAMCA(I) 

CONTINUE 

XLAMOA (1)31*000 

AT=SSC  (l)*OSC  AL*SQRT  (Y  (1,1)  ) 

IF  (IPR.GE.l)WRITc  (6,10  6)  (XLAMQA(I)  ,  1=1 ,  NM  AT  )  ,  AT 
NPP=N 

IF (ISIAS.NS.O )NPP=N-1 

CAll  3UILOR  (Y  , XLAMOA, NPP,(’AX) 

FORMAT  (5X,  *  SY  NT  ME  TIC  VECTOR,  ANO  SQRT  (Y  11  >*  ,/ ,  10 G 12 . 5 ) 

GO  TO  1000 

CONTINUE 

IF  ( ITPR, LE «  0) GO  TO  449 
WRITE (IUT, 44a ) 

CALL  PRTMAT  (Y  ,.N,N,MAX,0) 

00  442  I =1 , NM  AT 
DO  442  J=1,NMAT 
Z  ( I  ,  J)  =X  (I  +  1,  J't’l) 

CAcL  PROMAT  (Z  ,Y  ,NMAT,  10, MAX, 0  ) 

FORMAT  (2X,  9G1 1.  4) 

F0RMAT(2X,*  INVERSE  AND  PRODUCT  MATRICES*)  • 

CONTINUE 
XLAMOA (l)sl.O 
30  450  1=2, N 
XLAMOA (11*0. C 
30  450  J*i , NM AT 

XLAMOA ( I ) =  XLA MO  A  ( I ) -Y  ( I- 1  ,  J  )*  X  (J+l , 1) 

CONTINUE 

RETURN 

ENO 


SUBROUTINE 


MEQUAT 


PURPOSE: 


To  see  an  mxn  dimensional  matrix  B  equal  to  a  matrix  A 
of  the  same  dimensionality 


Equation:  B  *  A 


ROUTINE  VARIABLES:  M 
N 
B 
A 

NDIM 

IOPT 


Row  dimensionality  of  B  and  A 
Column  dimensionaltly  of  B  and  A 
Matrix  to  be  set 

Matrix  to  which  the  matrix  B  is  set 

Maximum  number  of  rows  permissible 

Print  option;  0  for  no  printing,  2  or  greater  for 
printing 


SU3 ROUTINE  MEQUAT  <M,N,S,A  ,NOIK,IOPn 


I0PT*3  SET  3  TC  ZERO 

1  3  EQUAL  TC  A 

IQ  3  TO  IOtNT  IT  Y 

DIMENSION  A  (NCIM,  1),3CNDIM,1> 
DO  33  1*1* M 
00  33  J=1.N 

IF (ICPT.Ni.l) 3(1* J)*0.0 


IF  (IOPT.t  J.iO  .AISO.I.EQ.J)  E  (I, J>  =1.0 
IF (I CRT. EQ.i) =(I»J)=A(I,J» 

CONTINUE 
R  E  T  v,  R  N 

end 


SUBROUTINE: 


PLOP 


PURPOSE:  To  plot  a  pair  of  columns  of  the  array  X 

(This  routine  may  be  substituted  by  user's  own  routine) 

SUBROUTINE  PLGP(NPT,NF,Y  ,NDIM  ,TQ,OT  ,LA3£l  ,IN0£P,  I5UF) 


KPTSNUMB  OF  TIME  PTS  (WARNING*  NOIM  SHOULD  0L.Gt.NPT  *2) 
NF=NLM3£R  OF  FUNS 

Y  (< ♦  )  DATA  ARRAY  OF  DIMENSION  NOIM.NF 
TOssINITIAi.  TIME,  0T=TIM£  INCREMENT 
LABEL,  INOEP  =  TITlES  FOR  Y  ANO  X  AXES 
DIMENSION  Y  (NDIM,  NF)  ,  YY(2)  ,lA6EL(1)  ,INDEP  <1 > 

DIMENSION  X  (512) ,I8UF(5i2) 

COMMON  /IO/IR.ILT, IPR, IT  PR , IZPR .IROUND, IPLT 

M*nF*NOIM 

Ml =M*1 

M2=M*2 

NPTl*NrT*l 

NPT2=NPT+2 

X (1)=T0 

□0  3  K=2,NPT 

X (K)=X (K-i) *0T 

DO  6  1=1, NF 

DO  d  K=NPT1,N0IM 

Y (< , I ) =Y (NPT, I) 

INITIALIZE  (LIQ.  INK,  12IN.FAPER)  MAX.  LE  NG  TM=60  IN 
CALL  PLOTMX  (6  G  .  0) 

SET  ORIGIN 

CAlL  PLOT (0  •«  -.  5,  3) 

CALL  FACTOR  (5  .0/6.5) 

BEGIN  PLOTTING 

CALi.  S C A Ll  ( X,  6 «  5  »  N  PT  »  1 ) 

CALL  SCALE (Y(l*l) ,10.0, M,l) 

CALL  AXIS(0«,  0.  .UPTIME  (SEC.), 

#-16»6»5»Q. ,X(NPTl) ,X  (NPT  2  )  ) 

CAL*  AXIS(0.»  0. , 1 6HRESPONSE S  Y  ,Y, 
•16,10.,,90.,Y(M1),Y(M2) ) 

WRITE  (6,6)  X  (NPT  1)  ,  X  (NPT2  ) 

WRITE  (6,7)  Y(M1)  ,Y  (M2) 

FCRMAT  (1X,*TC  ,OIV  (6.5  0IV)*»M1X,F7.3)  ) 

FORMAT  (IX, *YQ  ,01V  (10  DlV) * ,4  (IX, F7  .3  ) ) 

00  10  1=1, NF 
Y (NPT1 ,1) =Y (Ml) 

Y  (NFT2,I)=Y(M2) 

IF (I.EG.l)CAcL  LINE (X,Y( 1,1), NPT, 1, 1-1,1) 

IF  (I.  £C.  2)  CALL  CA  SHlN  (  X,  Y  ( 1  ,2  )  ,NPT,  1) 

CONTI.NLE 

CAtt  plot (io. ,o.,-3) 

RETURN 

END 
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SUBROUTINE: 


POLCON 


PURPOSE:  To  combine  the  factors  of  a  polynomial  in  order  to  produce 

the  coefficients. 


SUBROUTINE  PQlCCN(C»R2»K,N) 


A  PGUYINOMIAl  construction  program  needed  fcr  2T0S 
DIMENSION  C  <1),R2  (1) 

COMPLEX  C»  R2»  COMP 
DIMENSION  DC( 2) 

EQUIVALENCE  <C0MP,0C> 

NPl  =  N-*-l 
00101=2, NP1 
R2  <I)  =  0.0000 
R2  ( 1 ) =1 . 0D0  0 
0  0  *♦  I  =  1 ,  N 


C0MP=C (I) 

IF  (I.EC.K.OR.  (DC(1). EQ. 0.000. AND. OC  (2). EQ. 0.000)  )  GO  TO  «♦ 

002 JJ=1, I 

J=I-JJ+1 

R2 ( J+l ) =R2 ( J+ 1>  *C ( I) ♦R2 (Jl 
R2(1)=R2(1)*C<I) 

CONTINUE 

RETURN 

END 


EQUATION:  (x  -  c.)(x  -  c.)  -  (x  -  c  ) 

i  z  n 

,  ,  ,  n-1  .  n 

-  r,  +  r«x  +  ...  +  r  x  +  r  ,.x 

12  n  n+1 


VARIABLES:  C  Vector  containing  the  roots  of  the  factors 

R2  Vector  returning  the  coefficients  of  the  polynomial 
K  Exclude  Kth  factor  if  K^O 
N  Number  of  roots  contained  in  C 
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SUBROUTINE: 


POLRT 


PURPOSE:  To  find  the  roots  of  a  polynomial 

EQUATIONS:  a,  +  a„x  +  . . . .  +  a  ^,xn 

1  2  n+1 

<x  -  -P2-jq2)  —  (x  -Pn-jqn) 


ROUTINE  VARIABLES 

XCOF  Coefficients  of  the  polynomial  (XCOF(l)»a^) 

COF  Work  vector 

M  Order  of  polynomial 

ROOTR  Real  parts  of  the  roots  are  returned  in  this  vector 

ROOTI  Imaginary  parts  of  the  roots  are  returned  in  this  vector 

IER  Type  of  error,  if  any,  returned  in  this  integer  variable 


FURTHER  DESCRIPTION: 


SLaRCuTlN-i  FOLRT  (XCOF, CO F,M,P.COTR,  ROOT  I , IER ) 


COMPUTcS  THE  REal  ANO  CCMPLEX  ROOTS  OF  A  REAl  POlYnCMIAL 
DESCRIPTION  cf  parameters 

XCCF  -VECTOR  OF  M-l  COEFFICIENTS  OF  TbE  POLYNOMIAL 
ORDERED  FROM  SMALlEST  TO  lARGEST  POWtR 
CDF  -WORKING  VECTOR  OF  LENGTH  M+l 
M  -ORDER  OF  POLYNOMIAL 

ROOTR-RESULTANT  VECTOR  OF  LENGTH  M  CONTAINING  REAL  ROOTS 
OF  THE  POLYNOMIAL 

ROOTI-RcSUu TANT  VECTOR  OF  LENGTH  M  containing  THE 

CORRESPONDING  IMAGINARY  ROOTS  OF  THE  POLYNOMIAL 
IE  R  -ERROR  CODE  WHERE 
IER=0  NO  ERROR 

i er=i  m  less  than  cne 

I£R=2  M  GREATER  THAN  36 

Is.R=3  UNAoLE  TO  DETERMINE  ROOT  WITH  500  INTE  RATIONS 
ON  5  STARTING  VALUES 
IER  =  4  HIGH  ORDER  COEFFICIENT  IS  ZERO 

DIMENSION  XCO  F(i)  ,COF(l)  .ROCTR (1)  ,ROOTI  (1 ) 


OOU3LE  PRECIS  ION  >0,  YO  ,  X ,  Y ,  XP  fi  ,  YPR  ,  U  X,  U  Y ,  V ,  YT  , XT  , U  .  XT  2,  YT  2.SUMSQ, 

1  DX,OY, TEMP .ALPHA, XCO F ,C0 F , ROOTfi,RO OT I, ER 1, ER2, XSS , XS , Y SS  , YS , T Cl 
COMMON  /IO/IR,IlT,IFR,ITPR, IZPRtIROUNDf IFlT 

LIMITEO  TO  36TH  OROER  PCLYNOMIAL  OR  LESS. 

FLOATING  FOINT  OVERFLOW  MAY  OCCUR  FOR  HIGH  ORDEF 
POLYNOMIALS  9UT  WILL  NOT  AFFECT  THE  ACCURACY  OF  T HE  RESULTS. 


NEWTON-RAFHSON  ITERATIVE  TECHNIQUE.  THE  FINAL  ITtRATIONS 
ON  EACH  ROOT  ARE  PERFORMED  USING  THE  ORIGINAL  POLYNOMIAL 

rather  than  the  reouceo  polynomial  to  avoid  accumulated 

ERRORS  IN  THE  REOUCcD  PCLYNOMIAL. 


E  R  2  =  1 • 00+50 
TCu=l. 00-6 
IFIT=0 
N*M 
IER  =  Q 

IF  (XCOFCN+D)  10 .25,10 
10  IF(N)  15,15,32 

SET  ERROR  COCE  TO  1 


15  IE  R  =  1 

20  IF (IER)2Q0 , 20 1, 200 
WRITE (6,2C3)IER 

FORMAT (IX, *ERROR  CALLEO  FROM  FOlRT,  IER  =  *,I3) 
R  E  T  L  R  N 

SET  ERROR  COCE  TO  4 


25  I E  R  3 4 
GO  TO  20 

SET  ERROR  COCE  TO  2 


:2  (FCZF.F-l) 


30  IER  =  2 
GO  TO  20 

32  IF  (N-36)  35,3  5,  30 
35  NX=N 
NXX=Nfl 
N  2  =  1 

KJ1  =  N*l 
00  •* 0  L  =  i,KJl 
HT  =KJ1-l 

40  COF  <KT»=XCOF(L) 

set  initial  values 

45  XO=. 00500101 
Y0*0. 01000101 

ZtRC  INITIAL  VALUE  COUNTER 

- 5EGIN  ITERATION - 

X  ANO  Y  ARE  T hE  REAL  ANO  IMAG  PARTS  OF  ROOT 
I N  =  0 
50  X  =  XO 


INCREMENT  INITIAL  VALUES  ANO  COUNTER 


XO=-1Q.O*YO 

YO=-1Q.O*X 


SET  X  ANO  Y  TO  CURRENT  VALUE 

X  =  XC 

Y  =  YC 
IN=IN+1 
GO  TO  55 

55  I F I T=1 
XPR  =  X 
YPR  =  Y 

EVALUATE  POLYNOMIAL  ANC  DERIVATIVES 

59  I C  T  =  0 

60  UX=0.0 
UY=0.0 

Y  =0.0 
YT=0.0 
X  T  =  1 . 0 
U=COF(N»l) 

IF (L)  65,130, 65 

65  00  70  1=1, N 
U  =N-I+1 
TE  Mr*C0F („ ) 

XT  2  =  X*XT-Y*YT 

Y  T  2  =X*  YT  ♦Y  *  XT 
U=U>TEMP*XT2 
Vs V*T£MP*YT  2 
F  1  =  1 

UX»lX*FI*XT*TE'1P 
UY=LY-FI*YT*TE'1P 
X  T  =  X  T2 


:3  (FCLH T-2J 


70  YT=YT2 

SL'MSC  =  LX*JX*U  Y*UY 
IF(SLMSQ)  75,110,75 
75  OX=  (V*LY-U*UX  ) /SUt-SQ 
x=x+ox 

OY=-  (U«<JY-*V*UX)  /SLMSQ 
ILOC=75 

IF  ( ITPR.  GE  .  I)  WRITE (ILT, 442) IL OC 
YsY ♦DY 

xss=x 

YSS=Y 

IF  (YSS.EG.0.0DO)YSS  =  1.0O0 
IF  (XSS.EC. 0.0  DO )XSS=1. 030 
RMAGsSCRT(XSS*XSS*YSS*YSS> 

MOO IF. APR. 8  0  *SS  TO  * SS*. 0 00 Oi*RMA G  IN  NEXT  CARD 
£  R 1 2A3  S  (DX/(XSS+»  0000l*RMAG  ))  ♦  A  3S ( D  Y / < YSS  *.  0  0  001  ♦ft  MAG  ) ) 
IF (ITPR.GE. 1) WRITE  < ILT , 444 ) CX , XSS ,  0 Y , YSS , ER1.ER2 
IL0C=77 

IF  (ITPR.GE. 1)  WRITE  (ILT, 442)  IL  CC 

IF  (SRi.GT.ER2 )GO  TO  78 

ER2SER1 

XS=XSS 

YS=YSS 

IF  (ERI-TCl) 100, eG,8Q 

STEF  ITERATION  COUNTER 

80  ICT=ICT*1 

IF  (ICT-5  00)  60,85,85 
o5  IF  (IFIT)130,90,10  0 
MO  IF (IN- 5)  50 ,5  5, 55 

- cXIT  FROM  ITERATIONS - 


SET  ERROR  CODE  TO  3 

IER  =  3 
X  =  XS 

Y  =  YS 
ER1=£R2 

DO  105  L=1,NXX 

MT=KJl-w+l 

TEM?=XCOF(MT) 

XCOF(MT)  *C  OF  (  L  ) 

COF  (L) =TEM? 

ITEMFsn 
N  =  N  X 

N  X  =  ITEM? 

IF(IFIT)  120,  55,1  20 
IF(IFIT)  115,50,115 
XsXPR 

Y  s  Y  PR 
IFITsO 

I F  (A6S  (Y)-1.0D-8*A3S (X) >13  5,125,125 

A  L  P**  As  X  +  X 

SLMSQ=X*XfY*Y 

N  s  N  •  2 
GO  TO  1 4  C 

XsQ.Q 

NXsf.X-1 


,*,XX  =  NXX-1 


64  (PCLF.7-Z) 


135  y  =0 .0 

SUtiSQsQ.O 
ALPHA2  X 
N*N-1 

140  COF (2) »COF (2) ♦AuPhA*COF(l) 

1U5  00  150  L=2,N 

150  COFCi_>1IsCOF{l'*,1)  ♦  ALP  HA  «C  CF  (L  ) -SLMSQ*CO  F  C  L- 1) 

155  ROOT  I ( N2) 3 Y 
ROOTR (N2 ) =  X 

IF (ERl.GT.TCL ) WRITE (6, 554) N2.ER1 
554  FORMAT  (IX, ♦ERROR  ON  ♦  ,  13  ,♦  TH  ROOT  IS  *,010  ,3) 
£R2=1.00>50 
N2=N2*1 

IF(SLMSO)  160,165,160 
160  Y  =  -Y 

SUMSQ-0 • 0 
GO  TO  155 

165  IF(N)  20,20,45 
2  FORMAT  (2X, ♦TEST  IN  POLRT*,I5> 

4  FORMAT <1X,»0X ,X SS  ,  OY ,  YSS  , E R 1- £♦  ,2E  6 . 1 , 2 X, 2E8.1,2X, 2E5.2) 

RETURN 
ENO 


65  (POLF.T-4) 


SUBROUTINE: 

PURPOSE: 

EQUATIONS: 

ROUTINE  VARIABLES: 


PRDMAT 

This  subroutine  computes  the  product  of  two  square 
matrices. 

A  «-  A*B 


A,  B  N  x  N  matrices 
NDIM1  Maximum  row  dimension  of  A 

NDIM2  Maximum  row  dimension  of  B 


LOC 


An  integer  which  is  printed 


SUBROUTINE  FROM  AT  < A «  S  ,  N«  NCIf*l,NOIM2  ,LOC> 


COMPUTES  PROOUCT  A  =  A*o 

DIMENSION  A  (NOIM1  ,1>  ,3  (N0IM2,  1)  ,C  (10,10  ) 
IF  (wOC.GE.i)W  RITE  <6,51  LOC 
DO  31  1=1, N 
00  31  J  =  1 ,  N 
C (I,J) =0.0 
DO  21  <  = 1 » N 

C (I, J) =C (I, J)  +A (I ,K)*£  <K,J) 

CONTINUE 
DO  41  1=1, N 
DO  41  J  =  1 «  N 
A  (It  J)  =C  (I,  J) 

DO  45  1=1, N 

WRITE (6,1a)  (A (I ,J) , J*1,N ) 

F CRMAT  (*  L0CATICN/INTEGER=*,I5> 

FORMAT  (2X,  10G13.6) 

RETURN 

END 

FUNCTION  COM5  (N  ,M  ) 

CALCULATES  COM3INATION  M  CUT  OF  N 
IF(N.lE.O)GO  TO  99 
1  =  1 
L  0  =  1 

IF(M.EQ.0>GC  TO  6 

MN1=N-M*1 

DO  5  I=MN1,N 


L=L*I 

00  7  1=1,  M 
L  0=lD*  I 
C0M3=-/l0 
RETURN 
END 


SUBROUTINE:  prtmat,  prtvec,  prcvec,  prvec 

PURPOSE:  These  four  subroutines  perform  printing  of  arrays, 

PRTMAT  of  matrices  and  the  other  three  of  vectors. 
See  subroutines  for  comments. 

EQUATIONS : 


ROUTINE  VARIABLES 

Subroutine  PRTMAT 

A  Matrix  being  printed 
M  Its  row  dimensionality 

N  Its  column  dimensionality 

NDIM  Maximum  number  of  rows  permitted 
LOC  Use  LQOQ.  If  nonzeTO,  this  same  number  is 
printed. 


FURTHER  DESCRIPTION: 
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Si. 3  SOU  TINE  PR T MAT  ( A  , M  , N ,  NO  I M,  LOCI 


PRINTS  A  MATRIX,  AND  AN  INTEGER  (PERHAPS  A  LOCATION)  IF  LOC.GE.l 

DIMENSION  AtNOIH,  X» 

IF (LOC.GE. 1) WRITE (6,5) LOC 
DO  31  1=1, M 

WRITE (6, 13) (A (I,J) ,J*1,N) 

FGRMAT  (*  LOCATIGN/INTEGER**  ,15) 


FORMAT  (2X,  10G13  .6) 

RETURN 

END 


SUBROUTINE  PR  TVEC ( A, N) 


PRINTS  A  COMPLEX  VECTOR,  MAX  N*5 
WHEN  VARIABLE  IS  SINGLE  PRECISION 

COMPLEX  A (5 ) 

WRITE (6,1)  ( A( I) ,1  =  1, N) 

FORMAT  (1X,2E12.5,4(3X,2E12.5)  ) 

RETURN 

END 

SUBROUTINE  PRCVEC<A,N) 


PRINTS  A  COMPLEX  VECTOR,  MAX  N=5 
WHEN  VARIABLE  IS  OOU3LE  PRECISION 

COMPLEX  A 
DIMENSION  A  (1  ) 

WRITE (6,1)  ( A ( I) ,1*1, N) 

FORMAT (IX, 2012. 5,4 (3X, 2D 12.5) ) 

RETURN 

END 

SUBROUTINE  PRVEC(A,N) 


THIS  SUBROUTINE  OUTPUTS  DCL'SLE  PRECISION  SINGLE  DIMENSIONED  ARRAY 
DIMENSION  Ad) 

W S I T £  (  6 , 1 )  (  A  (  I)  ,1*1, N) 

FORMAT  (IX, 6016. 6) 

RETURN 

END 
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SUBROUTINE:  respon 


PURPOSE:  To  determine  the  response  of  the  digital  transfer 

function  H(z)  to  an  input  sequence  v(k).  The  coefficients 
of  H(z)  are  given  to  the  routine  in  the  NPNP2*»N+N+2  vector 
GAMMA 

EQUATIONS : 

be  *■  kf  *  '  *-  •  ■  •  +■  i,-.  i  ° 

Rw  *  - - - -o 

i  t-  o.f  a 1  f-  •  ■  •  l 


x 


K 


tc  v|,  t  •  *  ■  r  ^ 


ROUTINE  VARIABLES 


X  The  vector  which  returns  the  response  of  H(z) 

V  Vector  containing  the  input  sequence 

N  order  of  transfer  function  H 

GAMMA  Vector  of  coefficients  of  H 

GAMMA  *  (1,  a^ ,  •••  *  a^ )  ~b q i  "b ^ t  • • •  p  "b^) 

XLAMDA  Work  vector 

MP1  Number  of  response  points  generated 


FURTHER  DESCRIPTION: 


The  routine  assumes  zero  initial  conditions 
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SUBROUTINE  R£  SPON  (X,  V.  N,  GAMMA  ,  X..AM0  A  ,MPl> 


DIMENSION  X  ( 1  )  ,  V  (  1 )  ,  G  A  MM  A  ( 1 )  , XL AMD A (1) 

DOUBLE  PRECIS  ION  XSAV,AC,SO 

NMlsN-1 

NP1 =N+ 1 

NPNP1=N>N+1 

NPNP2=N*N+2 

00  19  1  =  1,  NPNP1 

XL  AMDA ( I ) =Q .0 

XS  A  V  =  G  #  0 

DO  20  K=i, MP1 

IF  (N.EC.  1)  GO  TO  25 

DO  21  1=1, NM1 

J=NP1-I 

XL AMDA (J>=XLAMOA(J-l) 

CONTINUE 
00  22  1=1, N 
J=NFNP2-I 

XL AMDA  CJ)=XLAMDA<J-1) 

XL  AMDA (1 ) =XSA V 
XLAMDA  (NPl  )  =V  «> 

XSAV=0 .0 

DO  23  I=1,NPNF1 

XSAV=XSAV“  GAMMA  <I*1)  ♦Xi.AMCA  (I> 

IF (0A3S(XSAV) .GE.1.0E10) XSAV=0.0 
X (Kl =XSAV 


SUBROUTINE: 


SIGNAL 


PURPOSE:  This  routine  generates  a  signal  which  is  a  weighted 

sum  of  exponential*sinusoid  terms 


EQUATIONS : 

m  -a  t 

f(t)  «  Z  W  e  1  Sin(B.t  +  <fi  ) 
i-1  1  11 


F(k)  -  f(kA) 


ROUTINE  VARIABLES 

P  Vector  returning  the  generated  signal 
NPT  Number  of  signal  points  generated 

AMP  Vector  of  weights  associated  with  each  exp*sinusoid  term 

SR  Vector  of  exponents  " 

SI  Vector  of  radian  frequencies  ” 

SPH  Vector  of  phases  " 

DT  Sampling  interval 
NCOMP  Number  of  terms 


FURTHER  DESCRIPTION: 

This  routine  is  useful  only  in  the  simulation 
mode  and  is  called  when  ISIM»2 
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SL'3  ROUTINE  SIGNAL  ( F,  NPT ,  A  MP  ,  S  R  ,  SI ,  S  PH ,  D  T ,  NC  OMP) 


DIrtENS  ION  F  (1  )  ,  AMP(l)  ,SR  <  1  >  ,S  I  ( 1)  ,  S  PH  ( 1 ) 

COMMON  /I3/IR,ILT,IPR, IT  PR, IZ PR .IROUND, IP  lT 

DOUBLE  PRECISION  A,3,C,X 

DO  12  K=1,NPT 

F  (0*0.0 

0  0  2  0  I  =  1,NC0MP 

A=SR(I)*DT 

a=SI(I)*OT 

C=SPH(I) 

00  15  K K s  1 ,  NP  T 
<=<K-1 

xsamp(I) 

IF  (A. Nc. 0.0 )X*X*OEXP(A*K) 

IF  (a.NE.O.Q)XsX*OSIN(3*K*CJ 
r (KK)=X*F(KK) 

CONTINUE 

IF (IPR.LT.2JG0  TO  30 
WRITE  ( ILT ,  9  I 

WRITE (ILT, 6) (F(<) ,K=1,NPT> 

WRIT£(I«-T,1) 

CONTINUE 
FORMAT (/> 

FORMAT (2 0 ( IX,  F5.2) > 

FORMAT (1 0 X , *  F  SIGNAL* > 

RETURN 

ENO 


SUBROUTINE: 


VEQUAT 


PURPOSE: 


To  equate  a  vector  Y  to  another  suitable  vector 


EQUATIONS :  Y (k) 


0  vector 
X(k) 

Y (k)*X(l) 

Y(k)+X(l) 

6k  (1,0, 0,0,..) 

1  (1,1, 1,1,..) 

Dimensionality  of  X 
The  vector  to  be  set 
Auxiliary  vector 
Not  used 

Option  parameter  (see  above) 


if  10PT-0 

-1  or  2  (print  also) 

-3 

-4 

-10 

-11 


VARIABLES: 


NPT 

Y 

X 

NPUL 

IOPT 


subroutine  vequat  <npt, y.x.npll, iopti 


lOPrsQ  Sc T  Y  TC  ZERO 

I  OR  2  SET  Y=X  (PRINT  IF  2) 

3  SET  Y=  Y*  CONST  X(l) 

4  Sc T  Y=  Y  +  CONST  X(l) 

9  SET  Y  TC  ZERO 

10  SET  Y= IMPULSE 

II  SET  YsSTEP 
DIMENSION  X  (1  )  »  Y  ( 1 ) 

IF  (I0PT.EQ.  0)  10 PT  =9 
00  33  K=1,NPT 

IF  ( IOPT. £3, 1, CR.IOPT.cQ.  2) Y (K )=X(K) 

IF  (ICPT.E J. 3) Y (K) *Y (K)*X (1» 

IF  ( ICPT.EQ.*) Y (K>  *Y(K) *X (1) 

IF  (IQPT.GE.9)  Y(K) =0.0 
IF  (I0PT. £‘3. 11  )Y  < K )  =1.0 
CONTINUE 

IF  (ICPT.EQ.2)  WRITE  (6,6)  (Y (X) ,K  =  1,NPT) 
FORMAT (EX, 10G 12.5) 

IF  (IOPT.E3.10  )  Y  (1>=1.0 

RETURN 

END 


SUBROUTINE: 


ZTOS 


PURPOSE: 


Given  the  z-domain  transfer  function  H(z)  this  routine 
computes  the  equivalent  s-domain  transfer  function  H(s) 


EQUATIONS: 


H(z)  H(s)  Impulse  invariant  if  IZTS=0 

Pulse  invariant  "  *1 

Trapezoid  invariant  "  »2 


ROUTINE  VARIABLES 

B  Vector  of  denominator  parameters 

A  Vector  of  numerator  parameters 

Comment  -  unfortunately,  these  names  are  opposite 
to  the  convention  used  in  the  rest  of  the  report- 
However,  the  reader  need  not  be  concerned  about 
this  unless  he  wishes  to  study  this  routine;  in 
the  latter  case,  he  should  bear  this  in  mind. 

N  Order  of  the  transfer  function 

DELTA  Sampling  interval 

IZTS  Option  for  type  of  conversion  (see  above) 


FURTHER  DESCRIPTION: 

See  Gold  and  Rader,  or  Oppenheim  and  Schafer,  or 
Stanley  for  the  theory  of  z-domain  to  s-domain 
converison. 
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SUBROUTINE  2TGS  Co  ,  A ,  N  *  D£  L  T  A  ,  I  ZTS) 


GIVEN  T  Hi  DISCRETE  DESCRIPTION  THIS  SUBRGL  TINE  COMPUTES  THE 
EGUIVAlENT  CONTINUOUS  DOMAIN  DESCRIPTION  OF  A  ulNEAR  0  YN  A  M I C 
system 

The  INPUT  ARRAYS  A  AND  3  ARE  FILLED  ACCORDING  TO  The.  DIFFERENCE 
EQUATION 

ECl  )*  Y  (<*  *5  (2>  *Y  U-l)  ♦.  ..  +  3  (N+1)*  Y(K-N> 

-A (1) *U (K)-A (2>*L(K-1) -...-A  (N*1)*U  (K-NI  =  0 
3(1)  MUST  tQUAL  1 

PClES  OF  THE  CONTINUOUS  ,'OMAIN  MUST  BE  DISTINCT  AND  NON-ZERO 
FOR  THE  TRNSFORMATICN  TO  Be.  VALID 


UPON  RETURNING  ARRAYS  A  AND  3  CONTAIN  THE  EQUIVALENT  COTINOODS 
DESCRIPTION  STCREO  ACCORDING  TO  THE  DIFFERlNTIAl  EQUATION 
6(1)*Y(TI+3C2) *OCl.Y(T) (N+l) *0 (N,  V(T)I 

-A(l>*U(T>-A(2)*D(l,U<T)>-...-A(N*l)*0<N,U(T))  =  0 
WhERE  0 ( M, F ( T )  )  =  THE  MTH  TIME  DERIVATIVE  OF  FUNCTION,  F 

3(N+i)  always  is  i 

N  =  ORDER  OF  SYSTEM 

N  (MAXIMUM)  =  ONE  LESS  THAN  THE  DIMENSION  SUBSCRIPT 


IZT  S=  0  -— >  IMPULSE 

IZTS=1 - -  PULSE 

IZTS=Z  - -  TRAPAZCIOAL 

delta  =  sampling  Interval  *  i/csampling  frequency) 

complex  CRtCA  ,CF,CG,CFl,CCNl,C0N2,C0NT 
COMPLEX  FOlECEQ ), ZRO (20) 

DIMENSION  3  (1  )  ,  A(l)  ,T£MP  (20  >  ,RR  <20  )  ,RI(2G  ),CR(20  >  ,CA(20  >  ,  CF(20>  , 
ICG (20)  ,CF1(20 ) • ZR  (20 )  ,  ZI  <  20  ) 

COMMON  ZIO/IR , ILT  » IPR, ITPR , IZPR t IROUND, IPLT 

CO NT-0,0 Du  0 

NF1=N*1 

IF  (A  (NP1I)  4  10  ,411,410 

ICHEC<=2 

GO  TO  401 

ICMEC<=1 

GO  TO  400 

CONT  =  A  (NF1)  ZB  (NPl) 

00*021  =  1, N 
A  ( I ) =A ( I ) -CON T*3  (  I) 

IF (ITPR. G£.l) WRITE (ILT, h41) 

CAlw  PClRT (5, TEMP,N,RR,RI,IER) 

IF (ITPR.GE.l) WRITE  (ILT  ,442) IE R 
FORMAT (2X, ’GOING  INTO  PCLRT*) 

FORMAT  (2X,*R£  TURNED  F R CM  POlPT,  IER=*,I2) 

3061*1,  N 

CR(I)=CMFlX(RR(I)  ,RI  (I)) 

CF  ( I)  *1.  030  0/CR  (I ) 

IF  ( IZPR.  G£.  1)  WRITE  (ft,  1002) 

F  C  R  M  i*  T  (*  ThE  POLE  S  OF  TnE  Z-OCMAIN*) 

IF (IZPR.GE.l) CALw  P RCVEC(CF.N) 

IF ( ;Z3R.C  0. 1) GO  TO  lie  1 
30  3  1  =  1  ,  '< 

3  C  N 1  =  1 . 0  >3  3  G  7o  (ZZGS~1) 


CON2=C.3DOO 
304 J  =  1 » N 

C0N2=CCN2*CR<  I)  ♦A  <N-J*i) 

IF<I-J)5,*,5 

C0N1=C0N1*  (1. 0000 -CR< I)*CF  ( J) ) 
CONTINUE 
CA ( I ) =CON2/CO  M 
IF <IZTS-1)225,224,230 
25  002221=1, N 

CA  ( I) =CA  (I)  /OtL  TA 
22  CR  (I)  =C*.OG  (CR  (I ) ) /OElTA 
GO  TO  226 
’4  0021=1, N 

CON1=ClOG  (CR<  1)  )/OELTA 
CA  (I)=CA  (I)*CR(I)*COM/(CR(I)  -l.ODOO) 
CR  (I>=COM 
GO  TO  226 


10  I Crt£CK=2 

002311=1, N 

C0N1=CL0G (CR( I> ) /DELTA 

CON  2=C  A  ( I)  *CR  (I)  /  <<1.000  Q-C  fid)  )*(  1.000  O-CR  (I))) 
C0NT=CCNT-C0N2*  <1.0000«-COM*0£LTA-CRCI)  ) 

CA  (I)=C0N2*C0M*CCN1#  OcL  T  A 
1  CR  (  I)  =C0M 
WRIT£<6,2131> 

31  FORMAT <2X,*TI M£  SCALED  FOR  SCATTERS R#) 

6  WRITE<6,1004)  IZTS 

0  4  FORMAT  (IX,  *S-  POlES  *,5X,-*SR*  ,9  X,«SI*  ,9X,*S  MAG*  ,9X,*FR*,7X, 
♦*IZTS=*, 12) 

OT=0. 001953125 
ROT  =DELT  A/OT 
R0T=1, 0 
DO  226  1=1, N 
SSR=-REAL(CR( I) )*ROT 
SSI  =  AIMAG  <CR( I)  )*  ROT 
SSM=CA9S  (CR  <1  ))  *RDT 
SSFR=(SSM/6. 26318  S3) 
e  WRIT£(6,1010) I,SSR,SSI,SSr,SSFR 
IF < IPR.Lt. 0 )G0  TO  1101 
10  FORMAT  (3X,I2,  5F12.4) 

IF  (IFR.SE.  DWRITE  (6,10  03) 

0  3  F  CRMAT  (*  NUMERATOR  CONSTANTS  OF  FACTORIZED  H(S)*) 

IF (IPR.GE.l)CALL  PRC VEC (CA , N) 

00  2 40  1=1, N 
1  POLE(I)=CR(I) 

C An.  PCi.CON(CR,CG»Q*N) 

0071=1, N?1 
CF  <I)=0.00  00 
009<=1 , N 

CAuw  R0lC0N(CR,CF1,K,N) 

009J  •  , N 

CF  (  J)  =CF  ( J  )  +C  F 1  <  J  )  *CA  <K) 

CF  (NP1  ) =0.000  G 
GO  TO  (403,40  4)  ,ICH£CK 
C  F (NP1 ) =C0N7 
CC405I=1 ,N 

CP<I)=C?C)4CONT*CG<I) 

CO  4  T  I  ,  Li  El 
r- (  ipr.le.i 


)GC  TC  520 


FIND  ZEROS 


00  507  1=1,  NP1 
07  A  ( I )  =CF  ( I ) 

nc  =  o 

0C  510  1=1, NFl 

IF  (A3S  (A  (I)  )  .GT.l.D-6)  NC  =  I-1 
1C  CONTINUE 

IF (NC.Ea.O) GO  TO  520 

N 1  =  N  0  ♦  1 
AK  =  A CN 1 ) 

00  515  1=1, N1 
L 5  A  (II  =A  (I )  /AK 

CALL  pclrt  (a,  temp,nc,rr,ri,ier) 

00  517  1=1, NO 
ZR  (I)=RR(I) 

ZI (I)  =  RI  (I) 

ZRO(I)=Cf*Pi.X(ZS  <I)  ,ZI(I) ) 

.7  CONTINUE 

IF (IFR.GE.2IWRITE (6,1007) AK 


C  7  FORMAT (*  ZEROS  OF  h(S),  NUMERATOR  CONSTANT 

IF (IPR,GE.2)CALL  P RT V£ C ( Z RO ,N 0) 

0  CONTINUE 

00201=1, NP1 
a (i)scG(ii 
A  ( 1 1  =CF  ( I) 

IF (IPR.Gt.l ) WRITE  (6,  10  05) 

05  FORMAT  (*  S-OOMAIN  DENOMINATOR*) 

IF  (IPR.3E.DCAu  L  PRVEC(3,NP1) 

IF (IPR.GE.l)WRITE (6,1006) 

C  6  F  0  RM  AT  (*  3-00  MAIN  NUMERATOR*) 

IF  ( IPR.GE.  1)CALL  PRVEC(A,NP1) 

01  RETURN 
ENO 


The  signal  considered  in  example  1  of  Section  IV  is 


x(k)  -  y(k)  +  w(k)  f 

-1  -2  -1  t 

where  y(k)  is  the  impulse  response  of  (1  -  1.92z  +  z  )/(l  -2.68z  +  , 

-2  -3  ! 

2.476z  -  0.782z  )  and  w(k)  is  additive  white  noise.  The  true  signal  ! 

% 

of  interest,  y(k),  is  shown  in  Fig.  B1  and  the  signal  under  test,  x(k),  is  ■ 

i 

shown  in  Fig.  B2.  ? 

-i 

Given  below  are  the  (card  deck)  input  to  program  POF-FILTER  and, 
succeeding* it,  the  printer  output  from  the  program.  ; 

INPUT  CARDS 


Fig.  B1  True  Impulse  response  of  a  third  order  transfer  function 
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APPENDIX  C 

MODELING  OF  A  SCATTERER  RESPONSE 

The  signal  of  interest  is  the  recorded  response  of  a  conducting 
pipe,  considered  in  Example  2  of  Section  IV.  On  site  digital  sampling  of 
the  response  was  carried  out  as  follows : 

k-0  to  k-300  Sampling  interval  «  0.390625  ns 

k«30l  to  302  "  «  0.703125  ns 

k-303  to  end  "  -  1.953125  ns 

For  analysis  purposes  we  resampled  the  first  300  data  points  by 
picking  up  every  5th  point;  the  remaining  data  were  used  as  such  to  gather 
a  total  of  245  data  points.  Note  that  we  have  ignored  the  intermediate 
sampling  rate  of  the  original  data  points  301-302. 

The  reconstituted  signal  was  shown  in  Fig.  12  by  the  solid  line. 

Two  runs  will  be  presented  below.  First  pertaining  to  the  signal 
obtained  above  (from  original  data).  The  second  pertains  to  differentiated 
(actually,  differenced)  signal  produced  from  the  resampled  response  (Fig.  13). 

Given  below  are  the  (card  deck)  input  to  program  POF-FILTER  and, 
suceeding  it,  the  printer  output  from  the  program;  first  for  the  signal 
itself  and  next  for  the  differentiated  signal. 

INPUT  CARDS  (RESAMPLED  SIGNAL) 
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PRINTER  OUTPUT  (conducting  pipe  response  analysis) 
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INPUT  CARDS  (DIFFERENTIATED  SIGNAL) 
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